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ABSTRACT 


The  development  of  a  universal  solution  of  the  main  problem  in  artificial  satellite 
theory  has  only  recently  been  accomplished  with  the  aid  of  high  powered  computers. 
The  solution  to  this  long  standing  problem  is  an  analytical  expression  that  is  similar  in 
form  to  the  two-body  solution.  An  analysis  is  presented  in  which  the  solution  is  com¬ 
pared  with  the  two-body  solution,  a  proven  numerical  solution,  and  actual  measured 
satellite  data.  The  solution  is  shown  to  be  significantly  more  accurate  than  the  two-body 
solution.  The  theoretical  accuracy  of  the  solution  is  confirmed.  The  solution  compares 
extremely  well  with  a  proven  numerical  solution  for  at  least  41  orbits  with  a  relative  error 
on  the  order  of  J20.  The  solution  compares  extremely  well  with  measured  satellite  data 
for  satellites  in  near  Earth  orbits.  For  a  satellite  in  orbit  at  an  altitude  of  approximately 
1000  kilometers,  the  solution  reduces  the  error  of  the  two-body  solution  by  about  95%. 
For  satellites  in  orbit  at  semisynchronous  and  geosynchronous  altitudes,  the  solution 
reduces  the  error  of  the  two-body  solution  by  at  least  50%.  The  solution  is  free  of 
singularities  and  is  valid  for  all  eccentricities  and  inclinations. 
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I.  INTRODUCTION 


With  the  development  of  any  new  method  or  theory  the  question  always  arises  on 
whether  the  approach  is  valid  or  practical  for  ordinary  use.  This  is  particularly  true  in 
the  prediction  of  satellite  motion.  Ever  since  Sir  Issac  Newton's  discovery  of  the  law 
of  universal  gravitation,  new  methods  have  been  developed  to  better  predict  the  motion 
of  the  heavenly  bodies.  Usually  the  method  contains  one  or  more  restrictions  that  limits 
the  practical  use  of  the  solution.  The  goal,  of  course,  is  to  develop  a  solution  that  is 
valid  and  possesses  no  restrictions.  Recently,  such  a  solution  has  been  formulated. 

This  analysis  continues  the  work  that  was  begun  by  Snider  [Ref.  l]  and  Sagovac 
[Ref.  2].  From  their  research,  a  higher  order  universal  solution  of  the  motion  of  an 
artificial  satellite  around  an  oblate  planet  was  developed.  The  solution  is  free  of 
singularities  and  is  theoretical  valid  for  all  orbital  parameters.  The  purpose  and  scope 
of  this  work  is  to  compare  the  solution  with  proven  numerical  solutions  and  actual 
measured  satellite  data  in  order  to  determine  whether  the  theoretical  work  is  valid  and 
practical. 

The  first  chapter  summarizes  the  development  of  the  theory  and  presents  the  sol¬ 
ution  in  its  entirety.  Also  included  is  a  somewhat  less  accurate  simplified  solution.  An 
explanation  of  the  solution  near  the  critical  inclinations  is  presented.  The  chapter  con¬ 
cludes  with  a  discussion  on  the  conservation  of  specific  mechanical  energy.  The  next 
chapter  describes  the  method  of  analysis  and  explains  the  type  of  integration  routine 
exercised  in  the  evaluation.  The  method  of  comparison  is  presented  next  in  which  the 
error  parameters  are  described  in  detail.  The  results  follow  which  include  both  a  detailed 
discussion  and  a  graphic  representation. 

The  analysis  is  objective  in  nature  and  designed  to  demonstrate  both  the  advantages 
and  disadvantages  of  the  theory  and  the  solution.  Before  the  solution  can  be  applied 
extensively,  a  general  understanding  of  its  strengths  and  weaknesses  must  be  determined. 
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II.  BACKGROUND 


A.  ORBITAL  KINEMATICS 

A  reference  system  for  a  planet  in  spherical  coordinates  (  r,  a,  /? )  is  shown  in  Figure 
1.  The  radial  distance  (r)  is  measured  from  the  center  of  the  planet  (O)  to  the  satellite 
(S).  The  line  (O  y)  is  in  a  direction  fixed  with  respect  to  an  inertial  coordinate  system. 
For  Earth,  the  line  (O  y)  is  in  the  direction  of  the  vernal  equinox.  The  right  ascension 
angle  (a)  is  measured  in  the  planet's  equatorial  plane  eastward  from  the  line  (0  y).  The 
declination  or  latitude  (/?)  is  the  angle  measured  northward  from  the  equator.  The  po¬ 
sition  vector  (r)  of  the  satellite  in  the  spherical  coordinate  system  is: 

r  =  r  ( cos  a  cos  /? )  bj  +  r  ( sin  a  cos  /? )  b2  +  r  ( sin  /? )  b3  (1) 

where  (  b, ,  b2 ,  h, )  are  orthonormal  base  vectors  fixed  in  the  direction  shown. 


polar  axis 


Figure  1.  Spherical  coordinate  system. 

A  reference  system  for  a  satellite  in  polar  coordinates  ( r ,  0)  within  a  rotating  orbital 
plane  is  shown  in  Figure  2.  The  satellite's  position  and  velocity  vectors  are  contained 
within  the  orbital  plane.  The  argument  of  latitude  (0)  is  measured  in  the  orbital  plane 


from  the  ascending  node  to  the  satellite.  The  inclination  (/)  of  the  orbital  plane  is  the 
angle  measured  between  the  equatorial  plane  and  the  orbital  plane.  The  longitude  of  the 
ascending  node  (Q.)  is  measured  from  the  line  (O  y)  to  the  ascending  node.  The  as¬ 
cending  node  lies  on  the  line  of  nodes  which  form  the  intersection  of  the  equatorial  and 
orbital  planes. 


orbi fal  plane 


Figure  2.  Orbital  plane. 

The  basis  (  b,  ,  b2 ,  h, )  may  be  transformed  into  another  orthonormal  basis 
(  B,  ,  B2 ,  Bj )  by  a  succession  of  three  rotations.  First  the  basis  (  b, ,  b2 ,  b, )  is  rotated 
about  the  b3  direction  by  the  angle  £1  The  basis  is  then  rotated  about  the  new  first 
coordinate  vector  by  the  angle  i.  The  final  rotation  is  about  the  new  third  coordinate 
vector  by  the  angle  0.  The  position  vector  (r)  has  only  one  component  in  the  rotating 
basis. 

r  =  rB,  (2) 
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The  components  of  r  in  the  fixed  basis  are: 

r  =  r ( cos  0  cos  Cl  —  sin  0  cos  /  sin  Cl )  b, 

+  r  ( cos  0  sin  Cl  +  sin-0  cos  /  cos  £2 )  b2  +  r  ( sin  0  sin  / )  b3 

Equating  the  components  of  equations  (1)  and  (3),  the  following  relations  among 
the  angles  (a  ,  /? )  of  the  spherical  coordinate  system  and  the  astronomical  angles  ( /, 
Cl ,  0  )  can  be  obtained. 


sin  /?  =  sin  0  sin  /  (4) 

cos  /?  =  cos  0  sec(a  -  £2)  (5) 

The  velocity  (dr/dt)  of  the  satellite  is  obtained  by  differentiating  equation  (2)  with 
respect  to  time  (r)  which  results  in: 

-£•  -  <« 

Equations  (2)  and  (6)  represent  the  orbital  kinematics  of  a  satellite  in  the  polar  co¬ 
ordinate  system.  The  position  and  velocity  vector  expressions  describe  the  motion  of  a 
satellite  in  a  particular  reference  system  and  provide  the  information  needed  to  develop 
the  equations  of  motion  in  that  system.  These  expressions  are  referenced  to  the  true, 
rather  than  mean,  orbital  plane  and  were  originally  formulated  by  Struble 
[Ref.  3,4,5]. 


B.  EQUATIONS  OF  MOTION 


The  motion  of  all  objects  is  mathematically  described  by  the  equations  of  motion 
that  govern  them.  For  an  oblate  planet,  the  expressions  for  the  kinetic  and  potential 
energies  per  unit  mass  of  an  orbiting  satellite  in  spherical  coordinates  are  respectively: 


V 


GM 

Y 


J,R 


1  +  (1  —  3  sin?7?) 


2  r 


(8) 


4 


In  the  above  equations,  (ftf)  is  the  mass  of  the  planet,  ( G )  is  the  universal 
gravitational  constant,  (R)  is  the  equatorial  radius  of  the  planet,  and  ( J2 )  is  the  coef¬ 
ficient  of  the  second  zonal  harmonic  of  the  planet's  gravitational  potential.  The  gov¬ 
erning  equations  of  motion  can  be  determined  by  substituting  equations  (7)  and  (8)  into 
Lagrange's  equations  which  are  represented  by: 


d_  8(T-  V) 

d‘  V  it) 


-4-(T-  V)  =  0 
oq 


(9) 


where:  q  =  a,  r,  and  /? 

Three  equations  result  from  Lagrange's  equations  which  describe  the  motion  of  the 
satellite.  The  three  equations  are: 


From  the  equations  of  motion,  two  integrals  result  which  are: 


r2cos2p~~-  =  constant  (13) 

T+  V  =  constant  (14) 

Equation  (13)  results  from  integrating  equation  (10)  and  equation  (14)  simply  states 
that  the  specific  mechanical  energy  of  the  satellite  remains  constant.  To  change  the  in¬ 
dependent  variable  from  t  to  0  ,  equations  (1) ,  (2) ,  and  (6)  are  used  in  conjunction 
with  equation  (13)  and  some  initial  conditions  to  form: 
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(15) 


dt_ 

dd 


r 2  cos i  ( 
,h0  cos  i0  V 


1  +  tan  0  cot  / 


di 

dd 


) 


Letting  u  =  p0ir  ,  and  using  equation  (15),  the  remaining  equations  of  motion  (11) 
-  (12)  can  be  rewritten  as: 


di 

d6 


2  Ju  sin  0  cos  0  sin  i  cos  / 


c. 


cos  / 


+  2  Ju  sin20  cos3/ 


(16) 


d2u 
dd 2 


u 


+ 


c2  cos2/ —  7c2//  cost 


2-  1*, 
;./  i 


—•  sin  6  cos  6  (3  cos2/  — 


472//2  sin30  cos  6  cos6/  (3  —  cos2/  )j-/{c4  +  47//c2  sin20  cos4/ 
4//272  sin40  cos8/ } 


7Q  _  tan#  di 
dd  sin  /  dt 


where: 


c  —  cos  i0 
s  =  sin/0 

7=  ^4 

2Po 


(18) 


Equation  (18)  results  from  uncoupling  the  equations  for  Q.  and  /.  The  angles  Q, 
and  /  can  be  uncoupled  by  applying  the  fact  that  the  orbital  p!ane  must  contain  the 
velocity  vector.  The  differential  equations  (16)  -  (18)  are  coupled  by  nonlinear  terms  and 
apparently  cannot  be  solved  analytically.  If  the  right  sides  of  equations  (16)  -  (18)  are 
expanded  in  a  Taylor  series  expansion  in  powers  of  J ,  the  equations  simplify  to: 


di_ 

dd 


•  *3 

—  27// sin  0  cos  0  sin /cos  /  ,  „  ,2  2  3  •  3«  rt  .  ,3x 

- 2 - +  47 // sc  sm  0  cos  0  +  0(7  )  (19) 
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COS2/'  .  JcOS2i  —  4«  sin2#  COS4/  .  2r,  .  _•  ,2rt  ,n  ..  .2. 

=  — r~H - ~ — < - 2 - • — f-  m  L I  -f-  sin  6  (7  cos  / 

c  c  |  c 

—  3)3  +  2u-~s'm  6  cos  0  (1  —  3  cos2/ )  —  2^-^-  sin20  cos2/ j 

+  ^4^(1-)^)- I]  PO) 

c  ’ 

+  3a  sin^-cos-/  +  A  sin  8  cos  9  [4-&A  +  sir2«  (j 

—  3  cos2/  )3+  ^  ~Hj‘  J2  sin20  cos2/ 1  +  (2(./3) 

~  =  -  2Juc  sin20 +  4J2u2c3  sin  0  +  O(J3)  (21) 


Each  of  the  neglected  terms  in  equations  (19)  -  (21)  are  indicated  by  the  (O)  symbols. 
These  are  terms  which  will  be  multiplied  by  J  to  the  third  power  or  higher.  Equations 
(19)  -  (21)  are  identical  to  those  used  as  a  stai  .ing  point  in  the  analysis  of  Eckstein,  Shi 
and  Kevorkian  [Ref.  63. 

C.  SOLUTION 

The  initial  breakthrough  of  an  analytical  solution  of  the  equations  of  motion  re¬ 
presented  by  equations  (15) ,  (19) ,  (20) ,  and  (21)  was  obtained  by  Danielson  and  Snider 
[Ref.  1,  73  .  Further  refinement  of  the  solution  w'as  later  formulated  by  Danielson, 
Sagovac,  and  Snider  [Ref.  2,  83  •  The  tlnee  authors  developed  the  solution  through  the 
extensive  use  of  an  algebraic  manipulation  computer  program  called  MACSYMA. 
Through  the  use  of  an  algorithm,  MACSYMA  was  able  to  solve  for  the  variables  u  ( 
or  r  ) ,  / ,  H  ,  and  dtldO  in  terms  of  0.  T'.c  solution  includes  all  terms  multiplied  by  J 
and  excludes  terms  of  order  J 2  and  higher.  In  : .  Yr  to  maintain  a  solution  of  order  J 
when  0  <  1/ J  ,  the  solutions  also  need  to  appropriately  include  terms  multiplied  by  J7d. 
The  solutions  which  analytically  demonstrate  relative  accuracy  of  order  J70  arc: 


7 


'  3s  ,  2. 
t - r—  +  e0  {. 


2 

e0 


(■~J4~)+i L(^-2)cos(2y-20) 


r  =  Pol\\  +  e0  cos  y  +  J 

<?0[15(2  +  <?0V  ~  W(4  +  <>oy  +  240  sinj  Jo(^~-2)  sin[0  +  coc] 

+  - - - - = - -  ~  - - - - - 

'  h.*'.<2-4) 

2 

+  -T*  E(6  +  5^q)s2  —  4(1  +  —  0q  +  (Oq)  4 — V-  cos (37g  —  ct)0) 


“  .  ^  i  r  2 

<?0V(15s2  -  14)  sin  70  (  -^r  •  • I  sinj  2c%  -  70 (  ~~  - 1) 


+ 


6(5  i' --4) 


e 2  e2 

+  -^l*  (9s2  -  4)  cos(y  +  30o  -  «o) +  ~f~  (6  ~  7i'2)  cos(y  +  &o~  wo) 


+  Y  C(5co  -  2)s2  -  2ejj ]  cos(>  +  0O  +  (o0)  +  e0s2  cos (0O  +  w0) 

2  2 

€(\  S  €ri 

-  cos(y  -  50o  +  3w0)  +  —  (3s  -  2)  cos(y  -  20o  v  2&>0) 

2  2  2 

e05  Cf\  -) 

cos(y  -  0O  +  3co0)  +  (35  -  2)  cos0;  -  30o  +  3o>0) 


16 

c0 


+  (I  —  3s2)  cos(y  —  20o)  +  ~-  (2  -  3s2)  cos  j’  +  s2  cos  20o 


(22) 


4 


(4  —  5s^)  cos(j.?  —  00  —  <w0)  +  ~y  (2s2  —  1)  cos(y  +  2  0O) 


16 

c0 


g2 

+  ~  (6  -  1  Is2)  cosO  +  20)  +  ~  (2  -  3s2)  cos(2>-  +  20) 


2 

+  ~~  [2^o  -  (2  +  5cq)520  cos  20  +  —■  (9s2  —  8}  cos  2 y 


3  e0s 


-  “~f~-  cos(y  ~  40o  J~  2<yo)  -  -f-  (-s2  +  1)  cos(v  +  2&>o) 

+  [2e02  -  (14  +  5,2)s2l  cosO  -  3  0O  +  w0)]  +  0(72,  J'O, ...)} 
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2 

y  =  e-ao  +  j{(-^— 2)(«-9o) 

4(1 12  -  75/  +  2 60s4  -  296s2)  *>[■«( -f~  i)  j  cos[V.n - 2) 


24(5s2  -  4)2 


+  /6 {  *°V(14  ~  15si)(1f  -  13?-C-°5-2^-  +  4  (15s2  -  13)  cos  26„  (23) 

1  24(5s2  -  4)  2 

2  2 

•r  ~~  (15s2-  13)  cos>{60  +  «o)  +  ~ -  (15s2  -  13)  cos(30o  -  ©0) 

+  “  [5(9eJ  +  34)s4  +  4^J  -  34)s2  -  56e2]}  +  0(j\  J20, ...) 

/  =  /0  +  sc7  -y  cos  2<?  +  "y  cost  ^  +  26)  +  -y  cos(y  —  20)  — y  cos  26q 

«?(14r-  15s2)sin[^(-^— 2)j  sin[^2<yo-70(^— 2) 

+  12(5s2  —  4) 

-  -y-  cos(30o  -  co0)  -  *y-  cos(0o  +  w0)|  +  0(J2,  J20, ...) 

Q.  =  £0  +  c7  jfl0  -  0  +  -y  sin  20  -  e0  sin  j;  +  -y-  sin(>  +  20)  -  -y  sin(>  -  20) 

eo(15s4  -  45s2  +  28)  sin|*./0(-y--  2)"]  cos|  2 w0  -  70(y; - 2) 

+  6(5s2  —  4)2 

-  -y  sin  2 00  +  e0  sin(0o  -  co0)  —  ■—  sin(3 0o  -  <o0)  ~  -y  sin(0o  +  <u0)}  (25) 

+  cJ20  {  —  ^  - —  cos  2<y0  -  <?0s2  cos(0o  +  G.'0j  +  -tj  (7s2  -  4) 

<-  12(5s2  —  4)  24 

-  — ——  cos(30o  ”  Wq)  _  -s2  cos  2^o  +  “12*  (6  —  •j2)}  +  ^(^2>  ’-) 


e2(14  r-  15s2)  sin 
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t  —  t0  + 


0  r  2  2 

>  +  J  /-2-|  1  +  J  —  y~  ~  cos  20  +  ^ofa2  -  1)  c°s y  +  -*jr-  cos(30o  -  ct>0) 

egs2(15s2 -  14) sin[y0(-^-  —  2) j  sin^2w0 -J9^~--2^ 


12(5s2  —  4) 


+  s2-l  (26) 


e0(l  -  2s  ) 


+  - -  cos(y  —  20)  +  c^s  20o  + 

+  ~~  COS (0O  +  £D0)  j  +  0(J2,  J2d,  ...)  jc# 


^o(3  ~  4s2) 


cos(y  +  2d) 


D.  SIMPLIFIED  SOLUTION 

As  shown  in  equations  (22)  -  (2c)  ,  if  6  is  restricted  such  that  6  <  1/7 ,  the  solution 
should  be  of  order  J  and  the  n  iglected  terms  should  be  of  order  J2 .  For  an  Earth  sat¬ 
ellite,  J  <  3/2  x  10"3 ,  so  for  at  least  100  revolutions  the  relative  error  should  be  on  the 
order  of  10"6 .  If  6  is  restricted  such  that  9  <  1  ,  all  of  the  terms  of  order  J20  in 
equations  (22)  -  (26)  can  be  neglected  if  a  relative  e  :re r  of  order  J  is  desired.  By  neg¬ 
lecting  terms  of  order  J70  ,  the  solution  simplifies  considerably  to: 

r  ~  ft/j1  +  c0  cosj^0  -  <u0  +  J(0-0O)(^ - 2) 

+  J  1  -  4-  e02(l  -  “j-  )  +  L2eJ  -  (2  +  5e02)s2]  cos  29 

2  2 

+  (9s2  -  4)  cos(0  +  3 0O  -  2c o0)  +  (6  -  7s2)  cos(0  +  90  -  2o>0) 

2 

+  -^T  (4  -  5s2)  cos {9  -  0O  -  2co0)  +  (2s2  -  1)  cos(0  +  2 90  -  (o0) 

2  2  e2 

-  -jg-  cos(0  -  0O  +  2co0)  +  (3s2  -  2)  cos(0  -  30o  +  2w0) 

2  2 

-  cos (9  -  500  +  2co0)  -f  —■  (3s2  -  2)  cos(0  -  2 0O  +  «0) 

+  ( 1  -  3s2)  cos(0  -  2 0O  -  cu0)  -f  — -  (2  -  3s2)  cos (0  -  w0) 


(27) 


+  -y  (9sl  “  8)  cos(20  -  2ft)0)  +  y-  (6  -  1 1  s2)  cos(30  -  w0) 

+  y  C(5^o  —  2)s2  —  2<?q]  co s(0  +  0O)  +  c0s2  cos(0o  —  w0) 

+  y*  (2  —  3s2)  cos(40  +  2coq)  +  y-  (3s2  —  2)  cos  2  <u0 

-  — g —  co$(0  -  46 0  +  (o0)-^-{s2+l)  co s{6  +  co0) 

+  -y  C2cq  —  (14  +  5cq)s2]  cos(0  —  30o)  +  s2  cos  2 0O 
+  \  [(6  +  5c02)s2  -  4(1  +  c02)]  co s(0  -  0O) 

2  -l 

+  — y—  cos(30o  —  <o0)  +  0{J2,  J20, ...)  j 
i  =  /0  +  scJ  j  cos  26  +  -y-  cos(30  -co0)  +  ~-  cos(0  +  co0)  -  —  cos  2 0O 

en  L  Cn  1  2  (28) 

-  -|-cos(3(?o-(wo)--|-cos(0o  +  «o)J  +  O(J2,/0I...) 

&  =  £20  4-  cJ  |0O  -  0  +  y  sin  20  -  e0  sin(0  -  <u0)  *f  ~  sin(30  -  co0) 

+  —  sin(0  +  co0)  -  —  sin  2 0O  +  e0  sin(0o  -  w0)  -  ~  sin(30o  -  co0)  (29) 
-  ^sin(0o  +  a>o)]  +  0(./2,./20,...) 


t  =  +  ~fj~  J  +  J  £ “ — 2 — “Cos 20  +  c0(s2  —  1)  cos(0  —  w0) 

c0(3  —  4s*)  en(  1  r-  2s2)  - 

+  - ^ - cos(30  -  a )0)  + - - - -  cos (0  +  a>0)  +  s2  -  1 

2  2 
c<)S  CoS 

+  —  cos(0o  +  (O0)  +  — y-  cos(30o  -  <y0) 


+  -y  cos  20o  -5-  <?(72, 730,  ...)p0 
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E.  THE  CRITICAL  INCLINATIONS 

As  shown  in  equations  (22)  -  (26)  ,  the  solution  appears  to  be  well  bounded  for  al¬ 
most  all  inclinations.  However,  two  particular  inclinations  immediately  appear  that  may 
produce  a  singularity.  A  possible  singularity  occurs  when  the  inclination  is  equal  to  ei¬ 
ther  63.43  degrees  or  116.57  degrees.  These  two  inclinations  have  produced  mountains 
of  literature  and  are  well  known  as  the  critical  inclinations.  However,  if  the  solution  is 
replaced  by  the  limit  as  the  inclination  approaches  either  critical  inclination,  the  solution 
remains  finite.  More  specifically,  the  solution  at  either  critical  inclination  is  as  follows: 

r  =  p0l\l+eQcosy  +  J  1  -^-  +  e02(l  + -|- (3s2  -  2)  cos(2v  -  20) 

+  ■—  [(6  +  5cq)s2  -  4(1  +  el)]  cos(>  -  0O  +  <w0)  +  ~ —  cos(30o  -  cu0) 

2  2 

+  (9 s2  -  4)  cosG'  +  30o  -  w0)  +  —-(6-  7s2)  cos(j/  +  0O  -  w0) 

+  y  C(5cq  -  2)s2  -  2el ]  cos(y  +  90  +  co0)  +  e0s2  cos (0O  +  co0) 

22 

e0S  €r\  i 

-  cos(y  -  50o  +  3 (o0)  +  —  (3 s2  —  2)  cos(y  —  20o  +  2co0) 

e2s2  2 

-  cos(y  -  0O  +  3co0)  +  — ■  (3s2  -  2)  cos(y  -  30o  +  3(u0) 

+  -y(l  -  3s2)  cos(y  -  20o)  +  -y(2-  3s2)  cosy  +  s2  cos  29 0  (31) 

2 

+  -j|-  (4  -  5s2)  cos(y  —  90  —  co0)  +  ~  (2s2  -  1)  cos(y  +  2 90) 

2 

+  yj"  (6—1  Is2)  CO  s(y  +  20)  +  YY  (2  -  3s2)  cos(2^  +  20) 

2 

+  "jy  [2 el  -  (2  -t-  5<?q)s2]  cos  20  +  -y|-  (9s2  -  8)  cos  2 y 
->  2 

OCnS  Cn 

-  — j—  cosO  -  40o  +  2<u0)  -  -y-  (s  +  1)  cos(y  +  2o)0) 

-r  ^  [2cq  -  (Id  +  5cq)s2]  cos(y- 30o  +  w0)J 
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[15(2  +  c02)s4  -  14(4  +  e0V  +  24]  sin(0  +  w0) 
yy  (15s2  -  14)  sin  2co0  +  0{J2,  J39, ...)} 

(32) 


(33) 

=  £>0  +  cJ  |0O  -  0  +  — ■  sin  26  -  e0  siny  +  y  sin(y  +  20)  —  sin(y  -  20) 

-  y  sin  2 0O  +  e0  sin(0o  -  w0)  -  y  sin(30o  -  co0)  -  —  sin (0O  -f  <y0)j  (34) 

2  2 
+  c720|  -jy  (6s2  -  7)  cos  2a)0  -  e0s2  cos (0O  -f  ty0)  +  y-  (7s2  -  4) 

2 

-  y-  cos(30o  -  «o)  -  s2  cos  20o  +  -y  (6  -  s2)}  +  0(J2,  J30, ...) 


2 

J-  =  +  - 2)(0 -<!„)]• 

2  2 
+  720{-y  (130s2  -  105s4  -  28)  cos  2  <u0  +  y-  (15s2  -  13)  cos  20o 

2  2 
€{\S  o  n 

+  — —  ( 1 5s  -  13)  cos(0o  +  <w0)  +  — y-  (15s2  —  13)  cos(30o  —  <y0) 
+  yr  [5(9cq  +  34)s4  +  4(9cq  -  34)s2  -  56c02]}  +  0(J2,  J39, ...) 

i  =  io  +  scJ  |  y  cos  20  +  y  cos (y  +  20)  +  —  cos (y  —  20) 

-  y  cos  2 0O  -  y  cos(30o  -  co0)  -  y-  COS(0O  +  (O0)y 

+  J2e  -y5 -  (14  -  15s2)  sin  2(0 0  +  0{J2,  J36, ...) 


-  'o+iH1+-'. 


(2 -3s2) 


2 - cos  20  +  c0(s  —  1)  cosj^  -r  s  —  1 


C0(I  ~  2s2) 


c0(3  —  4s2) 


2 - -  cos (y  -  20)  -J-  y-  cos  2 0O  -{-  ~^-y - cos(k  +  20)  (35) 


2  2  “I 

€qS  e0s  I 

+  ~~2~ co s(0o  4-  o>0)  +  — —  cos(3 90  -  co0) 

2  2 

+  J2d  (15s2  -  14)  sin  2 co0  4-  0(72,  J30,  ...)}*0 


Clearly,  equations  (31)  -  (35)  demonstrate  that  the  solution  is  indeed  finite  for  both 
critical  inclinations.  Equations  (31)  -  (35)  are  only  valid  for  the  critical  inclinations  and 
were  first  developed  by  Sagovac  [Ref.  2]  .  The  primary  purpose  in  developing  these 
equations  was  over  a  concern  in  computer  programming.  Some  computers  have  major 
problems  when  a  denominator  approaches  zero,  and  unlike  humans,  will  not  replace  a 
solution  with  its  limit.  Therefore,  depending  upon  the  accuracy  of  a  computer, 
equations  (22)  -  (26)  can  replace  equations  (31)  -  (35)  for  inclinations  near  the  critical 
inclinations.  It  should  be  noted,  how'evcr,  that  the  solution  itself  is  valid  and  bounded 
for  all  inclinations.  It  is  the  limitation  of  the  computer  that  creates  the  singularity. 

The  simplified  solution  w'hich  is  showm  in  equations  (27)  -  (30)  ,  is  valid  for  all  in¬ 
clinations.  Since  all  terms  of  order  Jjd  have  been  neglected,  the  troublesome  denomi¬ 
nators  mentioned  earlier  do  not  appear. 


F.  SPECIFIC  MECHANICAL  ENERGY 

For  all  satellites  under  the  influence  of  conservative  forces,  the  specific  mechanical 
energy  remains  constant.  Therefore,  ar.  ideal  analytical  check  of  the  solution  wrould  be 
to  see  if  indeed  the  specific  mechanical  energy  at  any  time  is  a  constant.  This  simple 
check  was  performed  by  Danielson,  Sagovac,  and  Snider  [Ref  I,  2,  7,  S3  by  substituting 
equations  (22)  -  (26)  into  equations  (7)  and  (S).  The  substitution  yields: 


r+K  _  GW  -  <g)  I  ~  3  sinVp)  .  ,  ^ 


2Po 


2  [W 


(36) 


The  first  two  terms  on  the  right  side  of  equation  (36)  represent  the  initial  specific 
mechanical  energy.  AH  other  terms  multiplied  by  J  in  equations  (22)  -  (26)  combine  to 
zero  when  substituted  into  equations  (7)  and  (S).  Equation  (36)  demonstrates  that  by 
neglecting  all  terms  of  order  J :  and  higher,  the  specific  mechanical  cncrgv  at  any  time 
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is  precisely  equal  to  the  initial  specific  mechanical  energy.  Obviously,  the  solution  sat¬ 
isfies  the  requirement  of  constant  specific  mechanical  energy. 
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III.  METHOD  OF  ANALYSIS 

A.  ORBITAL  PARAMETERS 
1.  Argument  of  Latitude  (  0  ) 

Figure  2  illustrates  that  the  position  of  a  satellite  at  a  particular  time  can  be 
described  by  the  argument  of  latitude  (0) ,  the  radius  magnitude  ( r )  ,  the  inclination  (/), 
and  the  longitude  of  the  ascending  node  (Q)  .  As  shown  in  both  the  solution  and  the 
simplified  solution,  r,  i,  and  Q  are  only  functions  of  J  and  the  argument  of  latitude  (0) . 
Since  J  is  a  constant  for  all  planets,  a  simple  determination  of  these  terms  is  trivial  once 
0  is  known.  However,  the  determination  of  0  is  not  trivial.  Although  it  would  be  ideal 
for  all  of  the  equations  to  be  analytical  expressions,  equations  (26)  ,  (30)  ,  and  (35) 
contain  an  integral  that  must  be  evaluated  in  order  for  0  to  be  determined.  Herein  lies 
the  key  to  the  solution.  Given  an  elapsed  time  between  observations,  how  can  d  be 
precisely  determined?  Since  the  initial  angular  momentum  ( h0 )  is  known,  this  term  can 
be  moved  to  the  left  side  of  equations  (26) ,  (30) ,  and  (35)  to  yield  equations  in  the  form 
of: 


(t  —  t0)/i0 


' o 

r\j,  0){1  +f(J,  6)}dd 


(37) 


If  r  was  not  a  function  of  0,  an  evaluation  of  the  right  side  of  equation  (37) 
could  easily  be  conducted  that  would  yield  an  analytic  expression.  However,  r  is  also  a 
function  of  0  and  the  only  practical  technique  in  evaluating  the  integral  is  through  nu¬ 
merical  means. 

Several  numerical  methods  could  be  used  to  evaluate  the  integral  depending  on 
the  speed  and  accuracy  one  desires.  Since  accuracy  and  not  speed  is  desired  in  this 
analysis,  a  Romberg  integration  routine  was  used  to  evaluate  the  integral.  Since  the 
right  side  of  equations  (26)  ,  (30)  ,  and  (35)  are  sinusoidal  in  nature,  the  Romberg 
•scheme  converged  quickly  and  accurately. 


16 


Since  9  defines  the  upper  limit  of  the  integral,  in  order  to  arrive  at  a  solution, 
an  initial  9  must  be  estimated.  Once  9  is  estimated,  the  integral  can  be  numerically  in¬ 
tegrated  and  the  result  can  be  compared  to  the  left  side  of  equation  (37).  If  the  com¬ 
parison  is  accurate  within  some  predetermined  error,  the  iteration  is  complete  and  6  has 
been  determined.  If  the  comparison  produces  an  error  that  is  unacceptable,  6  can  be 
incremented  either  up  or  down  and  the  integral  can  be  reevaluated.  Eventually,  the  it¬ 
eration  will  converge  and  6  will  be  determined.  An  algorithm  of  the  iteration  procedure 
is  as  follows: 

1.  Estimate  9. 

2.  Evaluate  the  integral. 

3.  Compare  the  result  with  the  left  side  of  equation  (37). 

4.  If  outside  the  limit,  go  to  (5).  If  within  the  limit,  stop. 

5.  Increment  9  up  or  down  as  needed,  go  to  (2). 

The  determination  of  9  involves  a  combination  of  two  errors.  The  first  error  is 
contained  in  the  numerical  evaluation  of  the  integral  itself,  while  the  second  error  in¬ 
volves  the  comparison  of  the  result  of  the  integration  with  the  left  side  of  equation  (37). 
Unfortunately,  the  errors  do  not  linearly  combine,  but  rather  multiply  since  the  numer¬ 
ical  evaluation  of  the  integral  is  inherently  nonlinear.  In  order  to  make  the  comparison 
error  meaningful,  the  evaluation  of  the  integral  must  be  made  as  precise  as  possible.  In 
order  to  avoid  determining  whether  an  error  is  due  to  computing  or  truncation  errors, 
the  numerical  technique  used  in  this  analysis  did  not  rely  on  a  step  size  constraint. 
Therefore,  the  relative  error,  in  general,  can  be  specifically  controlled.  Since  in  this 
analysis,  accuracy  and  not  speed  is  desired,  the  Romberg  integration  technique  was  uti¬ 
lized.  The  Romberg  technique  does  not  depend  on  any  specific  step  size  and  the  evalu¬ 
ation  of  the  integral  is  determined  through  a  converging  algorithm.  Also,  the  relative 
error  of  the  integration  can  be  specifically  controlled.  In  general,  the  relative  error 
normally  demanded  in  the  integral  evaluation  was  on  the  order  of  10“12,  and  the  relative 
error  of  the  comparison  was  on  the  order  of  10_I°.  Since  the  computer  program  utilized 
in  the  analysis  was  written  for  double  precision  accuracy,  these  types  of  relative  errors 
presented  no  significant  problems.  The  double  precision  accuracj  enabled  the  computer 
program  to  calculate  up  to  sixteen  digit  precision. 
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2.  Radius  Magnitude  ( r ) 


From  equations  (22) ,  (27) ,  and  (31)  ,  it  can  be  seen  that  the  radius  magnitude 
( r )  is  a  function  or  J  and  6.  Once  6  is  known,  r  can  be  evaluated.  From  the  appearance 
of  equations  (22)  ,  (27)  ,  and  (31) ,  it  is  not  obvious  how  r  will  behave  as  the  orbit  of  a 
satellite  progresses.  However,  from  observations  of  actual  satellite  motion,  it  is  clear 
that  the  orbit  should  behave  elliptically  with  r  varying  from  a  minimum  value  at 
periapsis  to  a  maximum  value  at  apoapsis.  The  magnitude  of  J  plays  an  important  role 
and  fortunately  for  most  planets,  oblateness  effects  act  as  a  pertuibation  in  comparison 
to  the  main  gravitational  force.  Therefore,  a  large  value  of  J  causes  larger  variations,  in 
r.  Since  equations  (22) ,  (27)  ,  and  (31)  contain  a  number  of  sine  and  cosine  terms,  a 
sinusoidal  behavior  should  be  expected. 

3.  Inclination  (  i ) 

The  solution  of  the  inclination  is  shown  in  equations  (24)  and  (33)  ,  and  the 
inclination  for  the  simplified  solution  is  shown  in  equation  (28).  In  general,  these  three 
solutions  are  quite  similar.  Again,  once  6  is  known,  i  can  be  evaluated  easily.  It  can 
be  seen  from  equations  (24)  ,  (28)  ,  and  (33)  ,  that  i  will  vary  slightly  from  an  initial 
inclination  as  the  orbit  of  a  satellite  progresses.  Also,  since  a  number  of  sine  and  cosine 
terms  are  present,  the  variation  should  be  sinusoidal  in  nature.  From  inspection  it  is 
clear  that  the  magnitude  of  the  variation  is  dependent  upon  the  magnitude  of  J  and  the 
initial  inclination.  The  variation  of  the  inclination  should  not  behave  in  a  diverging 
fashion,  but  rather  in  an  oscillatory  fashion  about  some  arbitrary  mean  inclination.  This 
behavior  is  consistent  with  observations  of  actual  measured  satellite  data.  The  driving 
factor  in  all  inclination  variations  is  the  magnitude  of  J.  Since  for  Earth,  J2  is  on  the 
order  of  10~3,  these  variations  should  be  quite  small. 

4.  Longitude  of  the  Ascending  Node  (  Q. ) 

The  solution  of  the  longitude  of  the  ascending  node  (Q)  is  shown  in  equations 
(25)  and  (34)  ,  and  the  longitude  of  the  ascending  node  for  the  simplified  solution  is 
shown  in  equation  (29).  As  expected,  all  solutions  are  quite  similar.  As  with  the  case 
of  r  and  /,  Q.  can  easily  be  determined  once  0  is  known.  Unlike  the  behavior  of  r  and 
/,  the  variation  of  Q.  is  very  predictable  and  highly  meaningful.  With  the  presence  of  9 
alone  in  equations  (25)  ,  (29) ,  and  (34)  ,  Q.  possesses  a  linear  relationship  with  0  and 
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as  6  increases  with  time,  the  variation  of  Q.  from  Q0  should  be  linear.  Depending  upon 
the  initial  inclination,  this  variation  will  be  either  positive  or  negative.  This  type  of  be¬ 
havior  is  clearly  consistent  with  the  classical  behavior  known  as  nodal  regression.  For 
an  oblate  planet,  nodal  regression  is  a  linear  property  whose  magnitude  and  direction 
depends  upon  the  radius  magnitude  and  inclination  of  the  satellite.  In  equations  (25)  , 
(29)  ,  and  (34)  ,  the  radius  magnitude  is  contained  in  the  J  term.  Therefore,  the  mag¬ 
nitude  of  the  nodal  regression  is  entirely  dependent  upon  the  magnitude  of  J.  From  the 
analysis  of  the  behavior  of  Q.  as  8  increases,  the  nodal  regression  behavior  should  be 
extremely  obvious. 


B.  ROMBERG  INTEGRATION  TECHNIQUE 


The  Romberg  integration  technique  is  a  powerful  integration  method  in  which  ar¬ 
bitrary  accuracy  can  be  achieved  in  a  relatively  efficient  manner.  The  method  combines 
any  type  of  relatively  inaccurate  quadrature  method  with  a  Richardson  extrapolation  in 
order  to  quickly  and  accurately  converge  on  a  solution.  In  this  analysis,  a  simple 
trapezoidal  quadrature  was  initially  used  to  estimate  the  integral  and  then  a  Richardson 
extrapolation  was  used  to  improve  the  integration  to  the  desired  accuracy  level.  The 
trapezoidal  quadrature  first  estimates  the  integral  with  a  single  interval.  The  estimate  is 
then  improved  by  using  2  intervals,  4  intervals,  8  intervals,  etc.  For  purposes  of  iden¬ 
tification,  the  results  can  be  labeled  /„' ,  /02 , 704 ,  and  so  on.  These  results  can  be  arranged 
in  column  form  in  preparation  for  a  Richardson  extrapolation  and  each  new  member 
represents  the  technique  of  halving  the  prior  interval.  The  length  of  the  column  is  de¬ 
termined  by  the  accuracy  that  one  desires.  Once  the  first  column  is  arranged,  a 
Richardson  extrapolation  can  be  performed  by  the  following  equation. 


/? 


tM/2 

h-  i 


The  values  of  /"can  be  arranged  in  tabular  form  as  shown  in  Table  1. 


(38) 
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Table  1.  A  SCHEMATIC  OF  ROMBERG  INTEGRATION 


/o 

- 

- 

Jo 

i 

- 

- 

- 

Jo 

it 

- 

- 

r8 

A) 

it 

il 

ll 

- 

716 

A> 

it6 

it 6 

IH 

z16 

To  test  For  convergence,  the  value  represented  by  if  is  compared  with  the  value 
represented  by  lf_x  .  If  these  two  values  are  within  some  predetermined  error,  then  if 
becomes  the  evaluation  of  the  integral.  IF  convergence  has  not  been  reached,  then  an¬ 
other  row  is  calculated  and  the  process  continues.  An  excellent  example  of  the  Romberg 
integration  technique  is  shown  in  Ferziger  [Ref.  9]  .  In  this  example  the  Following 
solution  of  the  integral  is  desired. 


I  = 


’  i 

exdx 

‘'o 


(39) 


From  elementary  calculus,  the  exact  solution  is: 

Iexact  =  2.718281828  (40) 

The  technique  of  Romberg  integration  of  the  integral  is  shown  in  Table  2.  The  rel¬ 
ative  error  of  /*  to  /2g  is  7.81  x  10-5 .  The  relative  error  of  /3  to  is  1.97  x  10"'° . 
As  can  easily  be  seen,  the  integration  is  converging  very  nicely  and  the  error  found  in  the 
Final  solution  is  less  than  the  error  demanded  within  the  Romberg  integration  scheme 
itself. 
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Table  2.  ROMBERG  INTEGRATION 


n 

*0 

Jn 

1 

1.859140914 

- 

- 

- 

2 

1.753931092 

1.718861152 

- 

- 

4  s 

1.727221905 

1.718318842 

1.718282688 

- 

8 

1.720518592 

1.718284155 

1.718281842 

1.718281829 

The  advantage  of  the  Romberg  integration  technique  over  a  simple  quadrature 
method  is  obvious.  The  number  of  intervals  that  must  be  evaluated  is  very  small  and 
the  relative  accuracy  is  very  high.  In  order  to  attain  the  accuracy  that  the  Romberg 
technique  delivers,  the  trapezoidal  method  would  need  to  divide  the  integral  into  several 
more  intervals.  This  would  be  highly  inefficient.  For  smooth  functions,  the  Romberg 
technique  is  very  effective  and  efficient.  Since  equation  (37)  is  sinusoidal  in  nature  and 
thus  relatively  smooth,  the  Romberg  integration  technique  was  used  to  evaluate  the  in¬ 
tegral.  If  equation  (37)  had  not  been  so  well  behaved,  another  integration  technique 
might  have  been  warranted.  The  Romberg  integration  scheme  is  the  heart  of  the  anal¬ 
ysis  and  can  be  found  in  the  computer  program  shown  in  Appendix  E. 
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IV.  METHOD  OF  COMPARISON 

A.  NUMERICAL  INTEGRATION  COMPARISON 

In  order  to  verify  that  the  theory  is  valid  for  practical  application,  the  solution  must 
be  compared  with  proven  numerical  solutions  and  measured  satellite  data.  By  compar¬ 
ing  the  solution  with  a  numerical  integration  of  the  equations  of  motion,  theoretical  ac¬ 
curacy  can  be  specifically  determined.  As  shown  in  equations  (22)  -  (26) ,  theoretically 
the  solution  is  accurate  to  order  J26  .  A  numerical  integration  comparison  will  determine 
whether  this  prediction  is  correct.  In  order  to  verify  the  solution,  the  following  param¬ 
eters  will  be  compared. 

1.  Delta  radius  vector  (|  Ar  | ) 

2.  Earth  arc  angle  (T) 

3.  Delta  omega  (AD) 

4.  Delta  inclination  (A i) 

5.  Delta  theta  (A0) 

6.  Delta  radius  relative  error  ( |  Ar  |  />•„) 

7.  Delta  theta  relative  error  (A0/0„) 

8.  Radial  track  error  ( RTE) 

9.  Along  track  error  {A  TE) 

10.  Cross  track  error  ( CTE) 

1.  Delta  Radius  Vector 

A  graphical  representation  of  the  delta  radius  vector  ( I  Ar  | )  and  the  Earth  arc 
length  (T)  is  shown  in  Figure  3.  The  delta  radius  vector  is  the  magnitude  of  the  vector 
separating  the  solution  radius  vector  (r)  from  the  numerical  integration  radius  vector 
(r„) .  Mathematically,  the  delta  radius  vector  can  be  expressed  as: 

Ar  =  r-rn  (41) 
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The  delta  radius  vector  describes  in  overall  terms  the  global  error  in  the  solution. 
Another  common  name  for  this  error  is  the  Euclidean  normed  difference  in  ephemerides. 
Although  the  delta  radius  vector  provides  ample  information  on  the  global  error  in  the 
solution,  this  error  can  also  be  expressed  by  a  different  parameter  that  will  be  called 
Earth  arc  angle. 


Figure  3.  Delta  radius  vector  and  Earth  arc  angle. 

2.  Earth  Arc  Angle 

Earth  arc  angle  (VF)  is  simply  the  angle  between  the  two  positions  if  viewed  from 
sea  level  on  Earth.  For  simplicity,  the  position  at  sea  level  was  chosen  such  that  the  arc 
angle  from  the  center  of  the  Earth  was  bisected.  By  using  the  law  of  sines  and  cosines, 
the  Earth  arc  angle  is  easily  determined.  Since  satellites  are  tracked  b>  instruments  on 
the  surface  of  the  Earth,  a  better  feel  for  the  global  error  can  be  attained  by  determining 
the  angle  between  the  two  positions.  Most  satellite  tracking  radars  possess  beamwidth 
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and  field  of  view  limitations;  therefore,  will  provide  useful  information  on  whether  the 
solution  is  accurate  enough  for  satellite  tracking  radars. 

3.  Delta  Omega,  Delta  Inclination,  Delta  Theta 

A  break  down  of  the  global  error  can  be  described  in  the  errors  of  delta  omega 
(AQ)  ,  delta  inclination  (A/)  ,  and  delta  theta  (A0)  .  As  mentioned  previously,  AQ  will 
provide  an  insight  on  the  motion  of  the  line  of  nodes,  and  specifically  if  nodal  regression 
is  present.  The  change  in  inclination  will  provide  information  on  the  movement  and 
stability  of  the  orbital  plane.  The  parameter  in  which  all  errors  are  based,  A 0  ,  will 
provide  much  information  on  the  source  of  the  global  error.  It  is  clear  that  small  errors 
in  Ad  will  contribute  significantly  to  the  accuracy  of  the  solution. 

4.  Relative  Errors 

The  verification  of  the  solution  will  lie  in  the  confirmation  of  the  relative  errors. 
The  delta  radius  relative  error  and  the  delta  theta  relative  error  will  demonstrate  the  ac¬ 
tual  accuracy  of  the  solution.  Both  parameters,  A r  and  AO  ,  demonstrate  a  theoretical 
error  of  J36  .  Therefore,  the  delta  radius  relative  error  should  be  on  the  order  of  J20, 
while  the  delta  theta  relative  error  should  be  on  the  order  of  J 2 .  Comparisons  of  the 
relative  errors  between  the  solution  and  the  numerical  integrati  an  solution  will  provide 
the  evidence  for  theoretical  confirmation. 

5.  Track  Errors 

Another  me  hod  to  break  down  the  global  error  is  in  terms  of  track  errors. 
Figure  4  shows  a  graphical  representation  of  radial  track  error  ( RTE)  ,  along  track  error 
(ATE) ,  and  cross  track  error  ( CTE) .  These  errors  can  better  be  described  by  referring 
also  to  Figures  1  and  2.  The  radial  track  error  is  the  error  in  the  radial  direction  or  in 
the  B,  direction.  The  along  track  error  is  the  arc  length  error  in  the  plane  defined  by  the 
solution  radius  vector  (r)  and  the  B.  direction.  Together,  these  errors  describe  errors  in 
three  orthogonal  directions  or  planes  as  compared  with  some  reference  position.  The 
reference  position  in  this  case  is  the  numerical  integration  solution.  A  mathematical 
derivation  of  these  errors  is  as  follows: 

Ar  =  r(r„  +  A r,  0  +  AO,  i  +  AO,  Q.  +  AQ)  —  f(>'n,  0,  i,  £2)  (42) 
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Using  equation  (i), 


.  Ar  =  (r„  +  Ar)(B1  +  AB,)-r„B1  (43) 

Ar  =  r„Bj  +  ArBj  +  ^ABi  +  ArABj  —  (44) 

Neglecting  higher  order  terms, 

Ar  =  ArBj  +  r„ABj  (45) 

Continuing,  defining  AB, , 

AB,  =  (B1.AB1)Bi  +  (B2.AB,)B2  +  (B3.AB,)B3  (46) 

Using  the  rotation  transformation  and  after  performing  considerable  algebra, 

B,.ABj  =  0  (47) 

B2.AB,  =  (A0  +  AQ  cos  f)  (48) 

B3  •  ABj  =  (A/  sin  0  —  AQ  cos  d  sin  i)  (49) 

Therefore,  using  equation  (45), 

Ar  =  (Ar)B,  +  r„(A0  +  AQ  cos  i)B2  +  r„(Aj  sin  0  —  AQ  cos  0  sin  /)B3  (50) 

From  equation  (50),  the  track  errors  can  easily  be  defined. 

RTE  =  r-rn  (51) 

ATE  =  r„(A0  +  AQ  cos  i)  (52) 

CTE  =  rn(Ai  sin  0  —  AQ  cos  0  sin  !)  (53) 
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A  graphical  representation  of  the  track  errors  is  shown  in  Figure  4. 


ATE 


Figure  4.  Track  errors. 

Examination  of  equations  (51)  -  (53)  demonstrates  that  the  track  errors  divide 
the  global  error  into  three  distinct  regions.  Radial  track  errors  obviously  describe  errors 
in  the  radial  direction.  Along  track  errors  are  similar  in  nature  to  Earth  arc  angle  errors, 
but  also  include  errors  due  to  nodal  regression.  Cross  track  errors  describe  orbital  plane 
errors  in  terms  of  both  inclination  errors  and  errors  due  to  nodal  regression. 

In  general,  all  these  parameters  should  give  an  excellent  insight  into  the  accu¬ 
racy  of  the  solution.  Also  included  in  the  numerical  comparison  will  be  the  simplified 
solution  and  the  two-body  solution.  The  simplified  solution  has  been  previously  pre¬ 
sented.  The  two-body  solution  can  easily  be  determined  by  simply  setting  7  =  0.  The 
analysis  of  the  numerical  comnarison  will  demonstrate  the  strensths  and  weaknesses  of 
all  the  different  solutions. 
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B.  MEASURED  DATA  COMPARISON 

In  order  to  observe  how  well  the  solution  models  actual  satellite  motion,  the  sol¬ 
ution  will  be  compared  with  actual  measured  data  from  operational  satellites.  To 
properly  evaluate  the  solution,  a  wide  range  of  orbit  characteristics  will  be  compared. 
These  characteristics  include  orbits  of  various  altitudes,  inclinations,  and  eccentricities. 
Also  included  in  this  comparison  will  be  the  simplified  solution,  the  two-body  solution, 
and  if  available,  some  particular  numerical  solutions. 

The  numerical  solutions  will  consist  of  two  forms.  The  first  is  a  numerical  solution 
that  only  includes  perturbations  involving  J2 ,  Jz ,  J4 ,  and  J5  .  The  second  numerical 
solution  will  include  the  following  perturbation  effects. 

1.  J2,J3,J4,  and  J5  . 

2.  Atmospheric  drag. 

3.  Sun  gravitational  effects. 

4.  Moon  gravitational  effects. 

5.  Solar  pressure  effects. 

From  the  analysis,  the  accuracy  of  the  solution  and  the  simplified  solution  can  be 
compared  to  a  numerical  solution  as  well  as  to  actual  measured  data.  The  weakness  of 
the  two-body  solution  will  also  be  demonstrated.  In  addition,  the  strengths  of  a  well 
modeled  numerical  solution  will  clearly  be  seen.  The  identical  error  parameters  described 
in  the  previous  section  will  also  be  used  in  the  measured  data  comparison.  From  this 
comparison,  the  advantages  and  disadvantages  of  the  solution  in  regard  to  actual  satel¬ 
lite  motion  will  clearly  be  demonstrated. 
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V.  RESULTS 


A.  NUMERICAL  INTEGRATION  COMPARISON 

To  verify  that  the  theory  is  valid  for  practical  applications,  the  solution  was  com¬ 
pared  with  a  proven  numerical  solution.  A  numerical  integration  computer  program 
called  UTOPIA  is  currently  in  use  at  the  Colorado  Center  for  Astrodynamic  Research 
(CCAR)  located  on  the  campus  of  the  University  of  Colorado,  Boulder,  Colorado. 
UTOPIA  is  primarily  used  to  model  a  wide  range  of  perturbations  and  can  predict  sat¬ 
ellite  motion  with  a  high  degree  of  accuracy.  The  UTOPIA  computer  program  was  de¬ 
veloped  at  the  University  of  Texas,  Austin,  Texas,  and  is  currently  in  use  at  several 
universities  and  research  centers.  The  solution  was  compared  with  the  UTOPIA  sol¬ 
ution  for  a  satellite  with  the  following  initial  conditions: 


rQ  =  7,3S6.IS  km 
iQ  =  90.03  degrees 
e0  =  0.003991 
0)Q  =  224.3S  degrees 
00  =  104.05  degrees 
Qq  =  322.63  degrees 
Itq  =  54,205. IS  km1  Is 
p0  =  7.371.29  km 
t0  =  0.00  seconds 

In  general,  these  initial  conditions  represent  a  slightly  retrograde  orbit  of  small  ec¬ 
centricity  at  an  altitude  of  approximately  1000  kilometers.  Essentially,  it  is  a  polar  orbit 
at  an  altitude  where  several  satellites  arc  currently  in  motion.  From  the  initial  condi¬ 
tions,  the  orbit  should  demonstrate  a  slight  easterly  nodal  regression.  But,  since  the  in¬ 
clination  is  so  close  to  90  degrees,  some  integration  routines  might  predict  zero  nodui 
regression.  In  this  comparison,  UTOPIA  only  modeled  the  J2  perturbation;  therefore, 
the  solution  should  compare  well  if  the  theory  is  valid.  All  error  parameters  depicted  in 
the  comparison  were  calculated  in  the  following  manner. 

A {ErrorParamctcr)  =  Theoretical  Solution  —  UTOPIA  Numerical  Solution  (54) 
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All  solutions  were  compared  at  one  hour  intervals  over  two  separate  periods  of  time. 
One  comparison  is  for  a  time  period  of  one  day  while  the  second  is  for  a  time  period  of 
three  days.  The  three  day  comparison  was  constructed  to  illustrate  the  effect  of  long 
term  errors  while  the  one  day  comparison  allows  for  a  more  detailed  analysis  during  the 
first  few  hours  of  motion.  The  period  of  rotation  for  the  satellite  is  about  105.26  minutes 
which  equates  to  approximately  13.68  orbits  per  day.  The  results  of  the  numerical  sol¬ 
ution  comparison  are  shown  in  Appendix  A. 

1.  Delta  Radius  Vector  Comparison 

The  comparison  of  the  delta  radius  vector  is  shown  in  Figures  5,  6,  and  7  in 
Appendix  A.  Figure  5  includes  the  comparison  of  the  solution,  the  simplified  solution, 
and  the  two-body  solution  to  the  numerical  solution.  If  the  solution  matched  the  nu¬ 
merical  solution  exactly,  the  delta  radius  vector  would  be  zero.  As  shown  in  Figure  5, 
the  solution  and  the  simplified  solution  match  extremely  well  with  the  numerical  solution 
while  the  two-body  solution  contains  gross  errors. 

A  more  detailed  plot  of  the  delta  radius  vector  comparison  is  shown  in  Figures 
6  and  7.  In  these  figures,  the  two-body  solution  is  excluded.  The  difference  in  the  sol¬ 
ution  and  the  simplified  solution  can  clearly  be  seen.  The  simplified  solution  produces 
a  diverging  sinusoidal  response  about  the  solution.  However,  up  to  approximately  four 
hours  of  motion,  the  solution  and  the  simplified  solution  are  nearly  identical.  The 
sinusoidal  behavior  of  the  simplified  solution  can  be  attributed  to  the  fact  that  all  terms 
multiplied  by  fO  have  been  neglected.  As  6  grows  with  time,  these  terms  become  sig¬ 
nificant  in  the  solution.  As  shown  specifically  in  Figure  7,  the  average  delta  radius  vec¬ 
tor  of  the  simplified  solution  clearly  diverges  from  the  solution. 

Figures  5,  6,  and  7  also  demonstrate  that  for  at  least  one  day,  the  delta  radius 
vector  for  the  solution  and  the  two-body  solution  are  nearly  linear  as  a  function  of  time. 
Other  comparisons  will  determine  whether  this  relationship  holds  true  and  will  be  shown 
later.  As  mentioned  earlier,  the  delta  radius  vector  is  a  global  error.  As  shown  in  Fig¬ 
ures  5,  6,  and  7,  the  solution  compares  well  globally  with  the  numerical  solution  and 
demonstrates  a  great  improvement  over  the  two-body  solution. 
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2.  Earth  Arc  Angle  Comparison 


The  comparison  of  the  Earth  arc  angle  is  shown  in  Figures  8,  9,  and  10  in  Ap¬ 
pendix  A.  From  inspection,  these  plots  are  nearly  identical  in  appearance  to  the  delta 
radius  vector  plots.  This  is  expected  since  both  the  delta  radius  vector  and  the  Earth  arc 
angle  represent  global  errors.  Figure  8  clearly  illustrates  the  large  error  generated  by  the 
two-body  solution.  The  two-body  solution  produces  unsatisfactory  long  term  satellite 
position  prediction.  After  one  day,  a  tracking  radar  would  have  a  difficult  time  detecting 
a  satellite  with  a  position  error  of  over  80  degrees. 

Figures  9  and  10  present  much  more  encouraging  results.  Again,  the  simplified 
solution  responds  in  a  sinusoidal  behavior  about  the  solution.  After  one  day,  the  posi¬ 
tion  error  of  the  solution  is  only  approximately  0.15  degrees.  Clearly,  the  solution  and 
the  simplified  solution  are  superior  to  the  two-body  solution.  Most  tracking  radars  can 
easily  handle  daily  position  errors  of  0.1 5  degrees.  In  general,  the  solution  and  the  sim¬ 
plified  solution  agree  very  well  with  the  numerical  solution. 

3.  Delta  Omega  Comparison 

The  comparison  of  the  delta  omega  angle  is  shown  in  Figures  11  and  12  in  Ap¬ 
pendix  A.  At  first  glance,  the  solution  and  the  simplified  solution  in  Figure  1 1  appear 
not  to  agree  well  with  the  numerical  solution.  However,  the  scale  of  delta  omega  is 
multiplied  by  10"4 .  The  numerical  solution  parallels  the  two-body  solution  nearly  ex¬ 
actly  and  predicts  almost  no  change  in  £2  .  In  other  words,  the  numerical  solution  pre¬ 
dicts  no  nodal  regression.  Easterly  nodal  regression  is  represented  by  a  positive  delta 
omega;  therefore,  it  is  clear  that  the  solution  and  the  simplified  solution  predict  greater 
nodal  regression  than  the  numerical  solution. 

On  a  larger  scale,  all  three  solutions  are  essentially  identical.  Since  the  initial 
inclination  is  so  close  to  90  degrees,  small  discrepancies  are  not  surprising.  The  delta 
omega  plot  does,  however,  invoke  confidence  in  the  solution.  Although  the  initial  in¬ 
clination  is  very  close  to  90  degrees,  the  solution,  the  simplified  solution,  and  the  nu¬ 
merical  solution  predict  easterly  nodal  regression.  This  result  is  significant.  Even  initial 
inclinations  close  to  90  degrees  produce  nodal  regression  in  the  correct  direction  for  the 
solution,  the  simplified  solution,  and  the  numerical  solution. 
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4.  Delta  Inclination  Comparison 


The  comparison  of  the  delta  inclination  angle  is  shown  in  Figures  13  and  14  in 
Appendix  A.  On  a  larger  scale,  all  of  the  solutions  compare  very  well.  On  the  scale 
shown  in  Figure  10,  the  two-body  solution  and  the  numerical  solution  are  nearly  iden¬ 
tical.  The  solution  and  the  simplified  solution  oscillate  about  an  error  of  approximately 
—2.0  x  10-5  degrees.  Obviously,  this  error  is  extremely  small.  In  general,  the  solution 
and  the  numerical  solution  agree  very  well. 

An  interesting  aspect  of  the  delta  inclination  comparison  is  the  sinusoidal  be¬ 
havior  of  the  solution  and  the  simplified  solution.  This  type  of  behavior  is  precisely  what 
was  predicted  in  the  earlier  analysis.  Figure  14  demonstrates  that  this  behavior  contin¬ 
ues  for  even  longer  periods  of  time.  On  a  larger  scale,  this  type  of  motion  would  not 
be  detectable. 

5.  Delta  Theta  Comparison 

The  comparison  of  the  delta  theta  angle  is  shown  in  Figures  15,  16,  and  17  in 
Appendix  A.  This  comparison  confirms  the  results  found  in  the  earlier  comparisons. 
The  two-body  solution  produces  very  large  errors,  while  the  solution  and  the  simplified 
solution  agree  very  well  with  the  numerical  solution.  Figures  16  and  17  again  illustrate 
the  typical  sinusoidal  response  of  the  simplified  solution  about  the  solution.  Since  the 
delta  theta  error  produces  all  other  errors,  the  excellent  results  found  in  the  earlier 
comparisons  are  now  not  surprising. 

6.  Delta  Theta  Relative  Error  Comparison 

The  comparison  of  the  delta  theta  relative  error  is  shown  in  Figures  18,  19,  and 
20  in  Appendix  A.  As  shown  in  Figure  18,  the  two-body  solution  demonstrates  a  rela¬ 
tive  error  of  2.37 ,  while  the  relative  errors  of  the  solution  and  the  simplified  solution  are 
much  smaller.  In  more  detail,  Figures  19  and  20  indicate  that  the  relative  error  of  the 
solution  is  2.S72.  This  result  confirms  the  theoretical  prediction  that  .ne  delta  theta 
relative  error  of  the  solution  would  be  on  the  order  of  J7.  Figures  19  and  20  illustrate 
that  initially  the  relative  error  of  the  simplified  solution  is  also  on  the  order  of  J2.  But 
as  6  increases  with  time,  the  relative  error  grows  in  a  sinusoidal  fashion.  This  result  is 
expected  since  the  simplified  solution  neglects  all  the  terms  multiplied  by  J70  .  In  gen- 
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eral,  the  results  shown  in  Figures  18, 19,  and  20  confirm  the  theoretical  relative  accuracy 
of  delta  theta  that  was  predicted  in  the  earlier  analysis. 

7.  Delta  Radius  Relative'  Error  Comparison 

The  comparison  of  the  delta  radius  relative  error  is  shown  imFigures  21, 22,  and 
23  in  Appendix  A.  As  shown  in  Figure  21,  the  two-body  solution  produces  a  relative 
error  that  isiinear  in  time  and  proportional  to  2.3  JQ .  Again,  the  relative  errors  of  the 
solution  and  the  simplified  solution  are  magnitudes  smaller. 

Figures  22  and  23  present  in  more  detail  the  relative  errors  of  the  solution  and 
the  simplified  solution.  The  relative  error  of  the  solution  is  very  linear  and  proportional 
to  2.8 J2d .  The  relative  error  of  the  simplified  solution  is  sinusoidal  in  nature  and  di¬ 
verges  from  the  solution.  However,  for  up  to  four  hours  of  motion,  the  relative  error 
of  the  solution  and  the  simplified  solution  are  nearly  indistinguishable.  Again,  the  re¬ 
sults  from  this  comparison  confirm  the  theoretical  prediction  that  the  delta  radius  rela¬ 
tive  error  of  the  solution  would  be  on  the  order  of  J2d . 

8.  Radial  Track  Error  Comparison 

The  comparison  of  the  radial  track  error  is  shown  in  Figures  24,  25,  and  26  in 
Appendix  A.  As  shown  in  Figure  24,  the  two-body  solution  oscillates  about  an  error 
of  approximately  -11.0  kilometers,  while  the  solution  and  the  simplified  solution  both 
produce  errors  that  are  dramatically,  smaller.  From  inspection,  the  two-body  solution 
also  appears  to  be  slowly  converging  as  time  increases. 

In  Figures  25  and  26,  the  solution  and  the  simplified  solution  produce  con¬ 
trasting  behaviors.  While  the  solution  remains  relatively  constant,  the  simplified  sol¬ 
ution  slowly  diverges  from  zero.  These  two  different  responses  continue  even  after  three 
days  of  motion.  Again,  the  neglected  J:9  terms  cause  the  significant  divergence  of  the 
simplified  solution.  Not  surprisingly,  the  solution  and  the  simplified  solution  are  clearly 
superior  to  the  two-body  solution. 

9.  Along  Track  Error  Comparison 

The  comparison  of  the  along  track  error  is  shown  in  Figures  27,  28,  and  19  in 
Appendix  A.  The  results  presented  in  this  comparison  parallel  the  results  found  in  the 
earlier  comparisons.  Since  the  inclination  is  so  close  to  90  degrees,  the  ACl  contribution 
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is  negligible  and  the  A 0  contribution  strongly  influences  the  responses.  As  a  result,  the 
along  track  error  comparison  is  practically  a  mirror  image  of  the  delta  theta  comparison. 

10.  Cross  Track  Error  Comparison 

The  comparison  of  the  cross  track  error  is  shown  in  Figures  30  and  31  in  Ap¬ 
pendix  A.  The  cross  track  error  is  strongly  influenced  by  A /  and  AQ .  Since  the  two- 
body  solution  produced  good  results  with  these  two  parameters,  it  is  not  surprising  that 
the  two-body  solution  agrees  well  with  the  numerical  solution.  Fortunately,  the  errors 
produced  by  the  solution  and  the  simplified  solution  are  also  very  small.  The  solution 
produces  a  maximum  cross  track  error  of  approximately  +  0.5  kilometers  after  one  day 
of  motion,  and  approximately  +  1.3  kilometers  after  three  days  of  motion. 

Clearly  in  this  comparison,  the  two-body  solution  is  superior.  However,  the 
large  errors  produced  by  the  two-body  solution  in  the  other  comparisons  easily  over¬ 
whelm  these  results.  In  global  terms,  the  two-body  solution  is  no  match  for  either  the 
solution  or  the  simplified  solution. 

B.  MEASURED  DATA  COMPARISON 

The  solution  was  compared  with  actual  measured  satellite  data  to  determine  the  al¬ 
titude  band  where  the  theory  works  best.  The  measured  satellite  data  was  obtained  from 
the  First  Satellite  Control  Squadron  (1SCS)  located  at  Falcon  Air  Force  Base,  Colorado. 
The  First  Satellite  Control  Squadron  tracks  several  satellites  for  the  Air  Force  and  was 
able  to  supply  measured  data  for  three  separate  satellites.  The  three  satellites  arc  cur¬ 
rently  in  motion  and  occupy  orbits  that  are  labeled  Near  Earth,  Semisynchronous,  and 
Geosynchronous,  respectively.  All  error  parameters  compared  in  the  earlier  numerical 
comparison  were  also  compared  in  this  comparison  using  the  measured  data  as  a  refer¬ 
ence. 

Included  in  all  the  comparisons  were  the  solution,  the  simplified  solution,  the  two- 
body  solution,  and  two  numerical  solutions.  The  two  numerical  solutions  were  also 
supplied  by  the  First  Satellite  Control  Squadion  and  are  labeled  Spacom  1  and  Spacom 
2,  respectively.  The  Spacom  1  solution  includes  all  perturbation  effects,  while  the 
Spacom  2  solution  only  includes  the  Earth's  harmonic  perturbations.  All  error  param¬ 
eters  in  this  comparison  were  calculated  in  the  following  manner. 
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A(Error Parameter)  =  Test  Solution  —  Measured  Data  Solution 


(55) 


Unfortunately,  the  First  Satellite  Control  Squadron  only  records  measured  data 
when  an  update  of  their  numerical  solution  is  required.  Routine  updates  are  usually 
conducted  after  about  seven  days  of  motion.  Therefore,  satellite  data  for  one  month 
usually  consists  of  only  four  data  points.  Although  more  data  points  are  needed  for  a 
more  detailed  analysis,  a  long  term  analysis  can  still  be  conducted.  The  analysis  of  each 
type  of  orbit  will  be  presented  separately. 

I.  Near  Earth  Orbit  Comparison 

The  near  Earth  orbit  comparison  possesses  the  following  initial  conditions. 


r0  =  7,776.58  km 
i0  —  98.81  degrees 
e0  =  0.0003071 
co0  =  9.57  degrees 
00  =  149.14  degrees 
=  37.10  degrees 

ho  =  53,664,37  km2js 
p0  =  7,224.89  km 
t0  =  0000Z  26  July  1990 

The  initial  conditions  of  this  satellite  represent  a  retrograde  orbit  of  small  ec¬ 
centricity  at  an  altitude  of  approximately  850  kilometers.  The  period  of  rotation  for  the 
satellite  is  about  101.89  minutes  which  equates  to  approximately  14.13  orbits  per  day. 
From  the  initial  conditions,  J2  should  be  the  dominant  perturbation.  The  orbit  should 
demonstrate  noticeable  easterly  nodal  regression.  If  the  theory  is  valid,  both  the  sol¬ 
ution  and  the  simplified  solution  should  agree  well  with  the  numerical  solutions  and  the 
measured  data. 

The  results  of  the  near  Earth  orbit  comparison  are  shown  in  Figures  32  -  43  in 
Appendix  B.  As  shown  in  the  figure the  solution  and  the  simplified  solution  agree  very 
well  with  both  the  Spacom  1  solution  and  the  measured  data.  The  fact  that  the  solution 
and  the  simplified  solution  produce  such  excellent  results  verifies  that  J2  is  the  dominant 
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perturbation  for  this  satellite.  Figures  32  -  43  also  demonstrate  the  larger  errors 
produced  by  the  two-body  solution.  In  almost  every7  comparison,  both  the  solution  and 
the  simplified  solution  are  far  superior  to  the  two-body  solution. 

One  surprising  result  is  the  poor  comparison  produced  by  the  Spacom  2  sol¬ 
ution.  In  every  comparison  the  Spacom  2  solution  either  models  the  two-body  solution 
exactly  or  produces  results  that  are  inferior  to  the  two-body  solution.  It  is  clear  that  the 
Spacorn  2  solution  does  not  model  the  Earth's  harmonic  forces  correctly.  An  explana¬ 
tion  for  the  poor  results  cannot  be  determined  in  this  analysis.  A  detailed  analysis  of  the 
force  modeling  in  the  Spacom  2  solution  must  be  completed  in  order  to  adequately  ex¬ 
plain  the  unsatisfactory  results. 

The  delta  omega  comparison  in  Figure  34  demonstrates  the  easterly  nodal  re¬ 
gression  produced  by  the  solution,  the  simplified  solution,  the  Spacom  1  solution,  and 
the  measured  data.  The  two-body  solution  represents  zero  nodal  regression.  Figure  35 
presents  the  delta  omega  comparison  at  a  much  smaller  scale  and  excludes  the  two-body 
solution.  In  this  figure,  much  more  detail  can  be  observed. 

There  is  only  one  comparison  in  which  the  results  are  mixed.  The  radial  track 
error  comparison  in  Figure  40  indicates  that  the  solution  produces  a  small  improvement 
over  the  two-body  solution  while  the  simplified  solution  actually  produces  a  greater  er¬ 
ror.  In  comparison  with  the  along  track  errors,  these  errors  are  small.  It  is  interesting, 
however,  that  the  radial  track  error  comparison  produces  such  mixed  results.  In  general, 
both  the  solution  and  the  simplified  solution  produce  results  that  are  in  excellent  agree¬ 
ment  with  the  measured  data  for  this  near  Earth  satellite. 
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2.  Semisynchronous  Orbit  Comparison 


The  semisynchronous  orbit  comparison  possesses  the  following  initial  condi¬ 
tions. 


r0  =  26,407.70  km 
i0  =  63.66  degrees 
e0  =  0.005860 
a)(j  =  318.19  degrees 
dQ  —  328.49  degrees 
£20  =  92.13  degrees 

Hq  =  102,892.59  km  Is 

Pq  =  26,559.96  km 

t0  =  0000Z  22  March  1990 

The  initial  conditions  of  this  satellite  represent  a  direct  orbit  of  small  eccentricity 
at  an  altitude  of  approximately  20,000  kilometers.  The  period  of  rotation  for  the  satellite 
is  about  717.96  minutes  which  equates  to  approximately  2.01  orbits  per  day.  An  im¬ 
portant  aspect  of  the  orbit  is  that  the  initial  inclination  is  very  close  to  the  critical  incli¬ 
nation  of  63.43  degrees.  Although  the  initial  inclination  is  not  exactly  that  of  the  critical 
inclination,  an  evaluation  of  the  solution  and  the  simplified  solution  near  this  important 
inclination  can  be  made.  From  the  initial  conditions,  the  orbit  should  demonstrate 
substantial  westerly  nodal  regression.  Also,  at  this  altitude,  the  dominance  of  the  J2 
perturbation  should  be  diminished.  Other  perturbations  that  are  not  modeled  should 
make  a  considerable  contribution  to  the  errors  in  the  comparison.  If  the  theory  is  valid, 
both  the  solution  and  the  simplified  solution  should  show  a  great  improvement  over  the 
two-body  solution. 

The  results  of  the  semisynchronous  orbit  comparison  are  shown  in  Figures  44  - 
53  in  Appendix  C.  As  predicted  earlier,  the  solution  and  simplified  solution  produce 
results  that  are  superior  to  the  results  produced  by  the  two-body  solution.  Figures  44 
and  45  present  the  global  errors  of  all  the  solutions.  In  global  terms,  the  solution  and 
the  simplified  solution  reduce  the  error  of  the  two-body  solution  by  nearly  one  half.  In 
effect,  the  J2  perturbation  accounts  for  approximately  one  half  the  error  produced  by  the 
two-body  solution.  The  remaining  error  which  is  represented  by  the  solution  and  the 
simplified  solution  is  caused  by  other  perturbing  forces.  Unfortunately,  the  results  of  the 
Spacom  2  solution  were  not  available. 


36 


The  delta  omega  comparison  in  Figure  46  demonstrates  the  easterly  nodal  re¬ 
gression  produced  by  the  solution,  the  simplified  solution,  the  Spacom  1  solution,  and 
the  measured  data.  Again,  the  two-body  solution  represents  zero  nodal  regression.  It 
is  clear  that  at  this  altitude,  the  J2  perturbation  produces  the  majority  of  the  nodal  re¬ 
gression.  The  delta  inclination  comparison  in  Figure  47  indicates  that  the  solution  and 
the  simplified  solution  produce  results  that  are  not  much  better  than  the  results 
produced  by  the  two-body  solution.  However,  the  error  after  30  days  of  motion  is  ex¬ 
tremely  small.  On  a  larger  scale,  the  solutions  would  seem  identical.  Since  the  inclina¬ 
tion  is  very  near  the  critical  inclination,  these  results  produce  more  evidence  in  support 
of  the  theory.  Clearly,  the  solution  and  the  simplified  solution  are  bounded  at  this  in¬ 
clination.  The  delta  theta  comparison  in  Figure  48  demonstrates  that  the  majority  of  the 
error  produced  by  the  solution  and  the  simplified  solution  originates  in  the  delta  theta 
error.  It  is  clear  that  the  two-body  solution  underestimates  the  value  of  0  while  the 
solution  and  the  simplified  solution  overestimate  the  value  of  0  . 

The  relative  error  comparisons  are  shown  in  Figures  49  and  50.  While  the  delta 
theta  relative  error  for  the  solution,  the  simplified  solution,  and  the  two-body  solution 
is  approximately  15.0  x  10"6,  the  relative  error  produced  by  the  Spacom  1  solution  is  far 
superior.  This  result  is  expected  since  the  Spacom  1  solution  models  several  more  in¬ 
fluential  perturbations.  The  delta  radius  relative  error  comparison  again  demonstrates 
in  global  terms  the  amount  of  improvement  that  the  solution  and  the  simplified  solution 
provide  over  that  of  the  two-body  solution. 

The  track  error  comparisons  in  Figures  51-53  produce  mixed  results.  While 
the  two-body  solution  produces  less  radial  and  along  track  errors,  the  solution  and  the 
simplified  solution  produce  much  less  cross  track  error.  In  comparison  with  the  along 
track  and  cross  track  errors,  the  radial  track  errors  are  small.  The  poor  results  produced 
by  the  solution  and  the  simplified  solution  in  the  along  track  error  comparison  is  due 
primarily  to  the  large  error  in  AO  .  The  very  large  error  produced  by  the  two-body  sol¬ 
ution  in  the  cross  track  error  comparison  is  due  primarily  to  the  very'  large  error  in  A£1 

In  summary,  although  the  solution  and  the  simplified  solution  arc  superior  to 
the  two-body  solution,  the  Spacom  1  solution  models  the  satellite  motion  more  precisely. 
However,  the  primary  reason  that  the  solution  and  the  simplified  solution  are  superior 
to  the  two-bodv  solution  is  due  exclusively  to  a  better  modeling  of  nodal  regression  or 
the  angle  £1  It  is  clear  that  the  solution  and  the  simplified  solution  model  the  J2  per¬ 
turbation  extremely  well.  The  Spacom  1  solution  is  expected  to  perform  better  since  it 
models  more  perturbing  forces. 
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3.  Geosynchronous  Orbit  Comparison 


The  geosynchronous  orbit  comparison  possesses  the  following  initial  conditions. 
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The  initial  conditions  of  this  satellite  represent  a  direct  orbit  of  small  eccentricity 
at  an  altitude  of  approximately  35,800  kilometers.  The  period.of  rotation  for  the  satellite 
is  about  1436.69  minutes  which  equates  to  approximately  1.00  orbit  per  day.  Since  the 
initial  inclination  is  slightly  greater  than  zero,  the  orbit  should  demonstrate  westerly 
nodal  regression.  However,  since  the  altitude  is  so  large,  other  perturbing  forces  that 
are  not  modeled  may  influence  nodal  regression.  At  a  geosynchronous  altitude,  the 
magnitude  of  other  perturbing  forces  approach  that  of  J2 .  Since  at  this  altitude  the  ef¬ 
fect  of  J2  is  so  diminished,  some  comparisons  of  the  solution,  the  simplified  solution,  and 
the  two-body  solution  may  be  nearly  identical.  As  a  result,  the  theory  may  not  be  any 
better  than  the  two-body  theory  for  satellites  in  a  geosynchronous  orbit. 

The  results  of  the  geosynchronous  orbit  comparison  are  shown  in  Figures  54  - 
63  in  Appendix  D.  The  global  error  comparisons  are  shown  Figures  54  and  55.  In 
global  terms,  the  solution  and  the  simplified  solution  produce  results  that  are  surpris¬ 
ingly  superior  to  the  results  produced  by  the  two-body  solution.  Evidently,  for  this  sat¬ 
ellite,  the  J2  perturbation  is  still  quite  dominant.  However,  the  other  comparisons  may 
present  a  different  picture.  Once  again,  the  Spacom  2  solution  generates  very  poor  re¬ 
sults. 

The  delta  omega  comparison  in  Figure  56  indicates  that  the  actual  perturbing 
forces  produce  easterly  nodal  regression.  Conversely,  the  solution  and  the  simplified 
solution  predict  westerly  nodal  regression.  It  is  obvious  that  other  perturbing  forces 
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influence  the  nodal  regression  of  this  satellite.  Although  the  Spacom  1  solution  is  su¬ 
perior,  even  this  accurate  numerical  solution  has  trouble  predicting  the  value  of  Q. .  The 
solution  and  the  simplified  solution  also  produce  poor  results  in  the  delta  inclination 
comparison  in  Figure  57.  All  solutions,  except  for  the  Spacom  1  solution,  produce 
identical  results.  Again,  on  a  larger  scale,  all  of  the  solutions  would  seem  nearly  identi¬ 
cal.  However,  this  detailed  analysis  does  demonstrate  a  weakness  in  the  theory.  The 
delta  theta  comparison  in  Figure  48  indicates  that  the  solution  and  the  simplified  sol¬ 
ution  are  inferior  to  all  solutions  including  the  two-body  solution.  Clearly,  other  per¬ 
turbing  forces  are  at  work. 

The  relative  error  comparisons  are  shown  in  Figures  59  and  60.  The  delta  theta 
relative  error  comparison  simply  repeats  the  results  found  in  the  delta  theta  comparison. 
However,  the  delta  radius  relative  error  comparison  is  much  more  reassuring.  Again,  in 
global  terms,  the  solution  and  the  simplified  solution  produce  better  results  than  the 
two-body  solution. 

The  track  error  comparisons  in  Figures  61  -  63  produce  mixed  results.  The  ra¬ 
dial  track  error  comparison  indicates  that  initially  the  Spacom  1  solution  is  inferior  to 
all  other  solutions.  However,  after  21  days  of  motion,  Spacom  1  is  the  superior  solution. 
Once  again,  the  radial  track  errors  are  small  when  compared  to  the  along  and  cross  track 
errors.  The  clue  to  the  favorable  global  results  of  the  solution  and  the  simplified  solution 
is  found  in  the  along  and  cross  track  error  comparisons.  The  solution  and  the  simplified 
solution  perform  much  better  than  the  two-bodv  solution  in  the  along  track  error  com¬ 
parison.  Although  the  two-body  solution  is  superior  to  the  solution  and  the  simplified 
solution  in  the  cross  track  error  comparison,  the  difference  is  small.  It  is  clear  that  the 
solution  and  the  simplified  solution  arc  superior  to  the  two-body  solution  due  to  a  much 
smaller  along  track  error. 

In  summary,  although  the  solution  and  the  simplified  solution  are  superior  to 
the  two-body  solution,  other  perturbing  forces  greatly  influence  the  satellite's  motion. 
At  this  altitude,  the  solution  and  the  simplified- solution  simply  do  not  model  the  satel¬ 
lite’s  motion  well.  Other  perturbing  forces  must  be  modeled  at  this  altitude  if  proper 
satellite  position  prediction  is  desired. 
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VI.  CONCLUSIONS  AND  RECOMMENDATIONS 


An  analysis  was  conducted  on  a  perturbation  solution  of  the  main  problem  in  arti¬ 
ficial  satellite  theory.  The  purpose  of  the  analysis  was  to  compare  the  solution  with 
proven  numerical  solutions  and  actual  measured  satellite  data  in  order  to  determine  if 
the  theoretical  work  is  valid  and  practical.  From  the  analysis,  the  following  conclusions 
can  be  made. 


1.  The  solution  and  the  simplified  solution  are  both  significantly  more  accurate  than 
the  two-body  solution.  The  relative  error  of  the  two-body  solution  is  on  the  order 
of  J6  while  the  relative  error  of  the  solution  and  the  simplified  solution  is  on  the 
order  of  J:9. 

2.  The  real  physical  effects  of  the  orbit  are  easily  distinguishable  in  both  the  solution 
and  the  simplified  solution. 

3.  The  solution  and  the  simplified  solution  compare  extremely  well  with  a  proven 
numerical  solution  for  at  least  41  revolutions  with  a  relative  error  on  the  order  of 
J26. 

4.  The  solution  and  the  simplified  solution  compare  extremely  well  with  actual  meas¬ 
ured  satellite  data  for  at  least  297  revolutions  at  altitudes  where  the  J2  perturbation 
dominates  ( e.g.,  near  Earth  orbits ).  For  a  satellite  in  orbit  at  an  altitude  of  around 
1000  kilometers,  the  solution  and  the  simplified  solution  reduce  the  error  of  the 
two-body  solution  by  approximately  95%. 

5.  The  solution  and  the  simplified  solution  compare  less  favorably  with  actual  meas¬ 
ured  satellite  data  at  semisvnchronous  and  geosynchronous  altitudes.  At  these  al¬ 
titudes,  however,  the  solution  and  the  simplified  solution  reduce  the  error  of  the 
two-body  solution  by  at  least  50%. 

6.  The  solution  and  the  simplified  solution  arc  free  of  singularities  and  are  valid  for 
all  orbital  parameters. 


Clearly,  the  solution  and  the  simplified  solution  model  the  J2  perturbation  very  well. 
The  equations  are  easy  to  implement  and  can  provide  quick  and  accurate  predictions  of 
satellite  motion.  However,  other  types  of  analytical  solutions  exist  that  are  more  accu¬ 
rate  than  the  solutions  described  here. 

One  such  solution  was  developed  by  Coficy  and  Alfriend  [Ref.  10]  through  re¬ 
search  that  was  conducted  by  Dcp.1.  CRef.  II 3,  Coffey  and  Dcprit  [Ref.  123,  and 
Alfriend  and  Coffey  [Ref.  133  •  The  solution  is  called  the  Analytic  Orbit  Prediction 
Program  generator  or  (AOPP).  Although  the  program  is  very  accurate,  AOPP  exten- 
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sively  utilizes  four  different  Hamiltonian  transformations.  As  a  result,  the  real  physical 
effects  of  the  orbit  are  not  easily  distinguishable. 

The  beauty  of  the  solution  and  the  simplified  solution  is  their  similarity  in  form  to 
the  well  known  two-body  solution  and  the  fact  that  a  satellite's  position  can  be  easily 
predicted  by  evaluating  only  one  integral.  Once  6  has  been  determined,  all  other  orbital 
parameters  can  be  calculated  easily.  The  structure  of  the  solution  and  the  simplified 
solution  is  ideal  for  implementation  with  onboard  spacecraft  computers. 

Before  the  solutions  can  be  adapted  for  practical  applications,  more  examination 
and  testing  of  the  theory  is  required.  In  order  to  provide  more  confidence  in  the  theory, 
the  following  recommendations  are  suggested. 

1.  The  solution  and  the  simplified  solution  need  to  be  compared  to  a  numerical  inte¬ 
gration  of  the  equations  of  motion  for  at  least  100  revolutions  to  confirm  the  the¬ 
oretical  accuracy  for  long  term  satellite  motion. 

2.  The  solution  and  the  simplified  solution  need  to  be  compared  to  several  more  di¬ 
verse  sets  of  actual  measured  satellite  data. 

3.  To  increase  precision,  the  solution  needs  to  include  the  higher  order  zonal  har¬ 
monics  of  the  gravitational  potential  (  e.g.,  Js ,  etc. ). 

4.  For  spacecraft  computer  implementation,  the  Lagrangian  coefficients  of  the  state 
transition  matrix  need  to  be  determined. 

For  onboard  spacecraft  navigation,  computers  make  use  of  the  state  transition  ma¬ 
trix.  Currently  the  Lagrangian  coefficients  of  the  two-bodv  solution  arc  the  only  matrix 
elements  that  have  been  determined.  An  excellent  formulation  of  the  two-body  state 
transition  matrix  is  shown  by  Battin  CRef.  143.  Once  the  Lagrangian  coefficients  of  the 
solution  are  developed,  onboard  spacecraft  navigation  can  be  greatly  improved. 
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APPENDIX  A. 


NUMERICAL  SOLUTION  COMPARISON  RESULTS 
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Figure  6.  .  Delta  radius  -vector  (1  day) 
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Figure  10.  Earth  arc  angle  (3  (lays) 
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DELTA  INCLINATION 
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Figure  18.  ,  Delta  theta  relative  error  (1  (lay) 
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Figure  19.  Delta  theta  relative  error  (1  day) 
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Figure  20.  ‘Delta  theta  relative  error  (3  days) 
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Figure  22.  .Delta  radius  relative  error  (1  day) 
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Figure  23.  Delta  radius  relative  error  (3  days) 
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Figure  24.  .  Radial  track  error  (I  day) 
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Figure  25.  Radial  track  error  (I  day) 
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Figure  27.  Along  track  error  (I  day) 
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Figure  29.  Along  track  error  (3  days) 
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Figure  3i.  Cross  (nick  error  (3  days) 
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Figure  33.  Earth  arc  angle  (21  days) 
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Figure  36.  Delta  inclination  (21  days) 
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Figure  37.  Delta  theta  (21  days) 
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Figure  38.  Delta  theta  relative  error  (21  days) 
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Figure  41.  ( Along  track  error  (21  days) 


TRACK  ERROR 


c 

c 

C 

o 

c 

o 

o 

U-j 

o 

U— ' 

-+-> 

13 

=5 

o 

00 

13 

o 

13 

O 

o 

CO 

Q 

Z 

LiJ 

O 

C 

00 

T) 

CO 

>N 

TJ 

o 

E 

CM 

E 

LlJ 

_J 

o 

1*-/ 

~Ol 

00 

o 

o 

o 

CJ 

C3 

E 

o 

o 

o 

O 

5 

CL 

CL 

00 

to 

F 

00 

CO 

1! 

II 

II 

II 

11 

© 

o 

<3 

<•> 

X 

m  310  WW3  >iOVdi  SSOdO 


,  HOURS 


TRACK  ERROR 


Figure  43.  ‘Cross  track  error  (21  days) 
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Figure  44.  *Dc!la  radius  vector  (30  days) 
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Figure  46,  •  Delta  omega  (30  days) 
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Figure  49.  Delta  theta  relative  error  (30  days) 
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Figure  50.  •  Delta  radius  relative  error  (30  days) 
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Figure  52.  .  Along  track  error  (30  days) 
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Figure  59.  Delta  theta  relative  error  (28  days) 
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I'iguive  61.  .Radial  track  error  (28  days) 
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Figure  63.  *  Cross  track  error  (28  days) 


APPENDIX  E.  COMPUTER  PROGRAM 
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PROGRAM  C0L02 
C 

C  ************************* 


c 

* 

* 

c 

* 

MAIN  PROGRAM 

* 

c 

* 

* 

c 

IMPLICIT  DOUBLE  PRECISION  (A-I,M-Z) 

CHARACTER*20  LINE 

DIMENSION  M( 100) ,MD( 100) ,U ’00) ,V( 100) ,WD( 100) ,0M( 100) ,OMD( 100) 
DIMENSION  I( 100) , ID( 100) ,F(  K'O) ,FD( 100) ,EC( 100) ,ECD( 100) ,A( 100) 
DIMENSION  R( 100) ,H( 100) ,N( 100) ,TH( 100) ,THD( 100) 

DIMENSION  RF( 100) , IF( 100) , IFD( 100) ,OMF( 100) ,OMFDC 100) ,THF( 100) 
DIMENSION  THFD( 100) ,P( 100) , JORBITC 100) ,DR( 100) ,DID( 100) ,DTHD(,100) 
DIMENSION  DOMD( 100) ,RX( 100) ,RY( 100) ,RZ( 100) ,RFX( 100) ,RFY( 100) 
DIMENSION  RFZ( 100) ,DRV( 100) ,ARC( 100) ,ARCD( 100) ,DAY( 100) ,HX( 100) 
DIMENSION  HY( 100) ,HZ( 100) ,VX( 100) ,VY( 100) ,VZ( 100) ,DT( 100) ,NX( 100) 
DIMENSION  NY( 100) ,NZC 100) ,RDVC 100) ,EXC 100) ,EYC 100) ,EZC 100) 
DIMENSION  NDEC 100) ,EDRC 100) ,VC 100) ,HTC 100) ,RDRFC 100) , INTAC 100) 
DIMENSION  ROMAC 100)  ,TH.’OC  100)  ,  ATE(  100)  ,CTEC  100) 

COMMON/OBLATE 1/DAY ,RX  RY,RZ ,VX,VY,VZ,DT,HX,HY,HZ,NX,NY,NZ,K,KK 
COMMON/ 0BLATE2/RDV ,  R ,  V ,  £X ,  P\ ,  EZ ,  MU.,  NDE ,  EDR ,  H ,  N ,  E ,  P ,  I ,  OM ,  W ,  F 
C0MM0N/0BLATE3/PI , EC , M , \ , HP , ER , TH , THD , RTD , MD , WD , OMD , ID , ECD 
C0MM0N/0BLATE4/FD , LINE , J , THF, THFD , IF , IFD , OMF , OMFD , RF , INT , ROM 
COMMON/OB LATE5/RJ , DR, Dll), DTHD , DOMD , ESTERR , ACTERR , TERROR , JVER 
C0MM0N/0BLATE6/RFX , RFY , RFZ , ARC , ARCD , RDRF , DRV , RJ2 , JN , JORBIT 
C0MM0N/0BLATE7/ INTA , ROMA , THJO , ATE , CTE 
C 

10  PRINTS 'ENTER  VERSION  C  SOLUTION  =  1,  SII-PLE  =  2,  TWO  BODY  =  3  )' 
READ*, JVER 

IFC  JVER.  EQ.  1.  OR.  JVER.  EQ.  2.  OR.  JVER.  EQ.  3) THEN 
GOTO  20 
ELSE 
GOTO  10 
ENDIF 
C 

20  PRINT*, 'ENTER  FIRST  POINT' 

READ*,K 

PRINT*, 'ENTER  FINAL  POINT' 

READ*,KK 

C 

PRINT1’1’,  'ENTER  RJ2' 

READ*,RJ2 

IFCRJ2.EQ.  1.  ODO)THEN 
RJ2=0.  00108263D0 
ENDIF 
C 

LINE=' . ' 

PI=3.  141592653589793238462643D0 
RTD=180. ODO/PI 
ER=6378. 137D0 
HTS=1.  0D0/806. 812D0 
MU=3.  986004 36D5 
C 

0PEN(UNIT=2.  STATUS='OLD' ,  FILE='/C0L02  OUT  A’) 
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0PEN(UNIT=3,  STATUS=’OLD’ ,  FILE='/C0L02  DSS  A') 
0PEN(UNIT=4,  STATUS=,OLD' ,  FILE=' /C0L020  DSS  B') 
C 

CALL  DATA 
CALL  ELEMENTS 
CALL  PINITIAL 
C 

J=1 

WRITE(6,6000)  ’POINT  =  ' , J 
WRITE(6,6000)  ’INTEGRATE  COMPLETED' 

WRITE(6,6000)  'FORMULA  COMPLETED’ 

WRITE(6,6000)  'INERTIAL  COMPLETED’ 

C 

DO  30  J  =  K,  KK 
WRITE(6,6000)  ’POINT  =  ’ , J 
CALL  INTEGRATE 

WRITE(6,6000)  ’INTEGRATE  COMPLETED' 

CALL  FORMULA 

WRITE(6S6000)  'FORMULA  COMPLETED’ 

CALL  INERTIAL 

WRITE(6,6000)  ’INERTIAL  COMPLETED’ 

30  CONTINUE 

CALL  RESULTS 
C 

CLOSE  (UNFD=2) 

CLOSE (UNIT=3) 

CL0SE(UNIT=4) 

C 

6000  FORMAT( 3X,A,I3) 

STOP 

END 


*:V********5'«V**5V****?V5V******A 


*  SUBROUTINE  DATA 

*  * 
*iViVV»V^5VynV:V:V*)V?V5'«'oV>Vi,«V>V5VVfi,r^)Vyr;V 


SUBROUTINE  DATA 

IMPLICIT  DOUBLE  PRECISION  (A-I,M-Z) 

CHARACTER*20  LINE 

DIMENSION  M( 100) ,MD( 100) ,E( 100) ,W( 100) ,WD( 100) ,0M( 100) ,OMD( 100) 
DIMENSION  I( 100) , ID( 100) ,F( 100) ,FD( 100) ,EC( 100) ,ECD( 100) ,A( 100) 
DIMENSION  R( 100) ,H( 100) ,N( 100) ,TH( 100) ,THD( 100) 

DIMENSION  RF( 100) , IF( 100) , IFD( 100) ,OMF( 100) , OMFD( 100) ,THF( 100) 
DIMENSION  T1IFDC 100) ,P( 100) , JORBITC 100) ,DR( 100) ,DID( 100) ,DTHD( 100) 
DIMENSION  DOMD( 100) ,RX( 100) ,RY( 100) ,RZ( 100) ,RFX( 100) ,RFY( 100) 
DIMENSION  RFZ( 100) ,DRV( 100) ,ARC( 100) ,ARCD( 100) ,DAY( 100) ,HX( 100) 
DIMENSION  HY( 100) ,HZ( 100) , VX( 100) , VY( 100) , VZ( 100) ,DT( 100) ,NX( 100) 
DIMENSION  NY( 100) ,NZC 100) ,RDVC 100) ,EXC 100) ,EYC 100) ,EZC 100) 
DIMENSION  NDEC 100) ,EDRC 100) , VC  100) ,HT( 100) ,RDRFC 100) , INTAC 100) 
DIMENSION  MONTHC 100) ,DATEC 100) ,KOURC 100) ,MINC 100) ,SECC 100) 
DIMENSION  ROMAC 100) ,THJOC 100) ,ATEC 100) ,CTEC 100) 

COMMON/OBLATE 1/DAY , RX , RY , RZ , VX , VY , VZ , DT , HX , HY , HZ , NX , NY , NZ , K , KK 
COMMON/OB  LATE  2 /RDV , R , V , EX , EY , EZ , MU , NDE , EDR , H , N , E , P , I , OM , W , F 
C0MM0N/0BLATE3/PI , EC , M , A , HT , ER , TH , THD , RTD , MD , WD , OMD , ID , ECD 
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COMMON/ 0BLATE4/FD , LINE , J,THF,THFD , IF , IFD ,OMF, OMFD ,RF , INT,ROM 
C0MM0N/0BLATE5/RJ , DR , DID , DTHD , DOMD , ESTERR , ACTERR , TERROR , JVER 
C0MM0N/0BLATE6/RFX , RFY , RFZ , ARC , ARCD , RDRF , DRV ,R J2 , JN , JORBIT 
COMMON/OBLATE 7 / I NTA , ROMA , TH JO , ATE , CTE 

READ  IN  EMPHERIS  DATA 

OPEN(UNIT=l,  STATUS=' OLD1 ,  FILE=*/C0L02  DAT’,  ACTION=’READ’ ) 
DO  10  J  =  1,  KK 

READ(1,*)  MONTH(J) ,DATE(  J) ,H0UR(J) ,MIN( J) ,SEC( J) 

READ(1,*)  RX( J) ,RY( J) ,RZ( J) 

RX( J)=RX( J)/1000.  0D0 
RY( J)=RY ( J) / 1000.  0D0 
RZ(J)=RZ(J)/1000.  0D0 
RE AD( 1,*)  VX( J) ,VY( J) ,VZ( J) 

VX(J)=VX(J)/1000.  0D0 
VY( J)=VY( J)/1000.  0D0 
VZ( J)=VZ( J)/1000.  0D0 
CONTINUE 

CONVERT  PARAMETERS 
DO  20  J  =  1,  KK 

DAY(J)=DATE(J)+((3600.0D0*HOUR(J)+ 

+  (60. 0D0*MIN( J)+SEC( J) ) ) ) / 86400.  0D0 

DT( J)=(DAY( J)-DAY( 1))*24.  0D0*3600.  ODO 
CONTINUE 

CL0SE(UNIT=1) 

RETURN 

END 


*5V*A*****>V5V**>V*?V****?Wn\nViVTiV?ViV*iV*5V 


*  SUBROUTINE  ELEMENTS 

* 


it 

J. 


SUBROUTINE  ELEMENTS 

IMPLICIT  DOUBLE  PRECISION  (A-I,M-Z) 

CHARACTER*20  LINE 

DIMENSION  M( 100) ,MD( 100) ,E( 100) ,W( 100) ,WD( 100) ,OM( 100) ,OMD( 100) 
DIMENSION  I( 100) , ID( 100) ,F( 100) ,FD( 100) ,EC( 100) ,ECD( 100) , A( 100) 
DIMENSION  R( 100) ,H( 100) ,N( 100) ,TH( 100) ,THD( 100) 

DIMENSION  RF( 100) , IF( 100) , IFD( 100) ,OMF( 100) ,0MFD( 100) ,THF( 100) 
DIMENSION  THFD( 100) ,P( 100) , JORBITC 100) ,DR( 100) ,DID( 100) ,DTHD( 100) 
DIMENSION  D0MD( 100) ,RX( 100) ,RY( 100) ,RZ( 100) ,RFX( 100) ,RFY( 100) 
DIMENSION  RFZ( 100) ,DRV( 100) ,ARC( 100) ,ARCD(100) ,DAY( 100) ,HX( 100) 
DIMENSION  HY( 100) ,HZ( 100) ,VX( 100) , VYC 100) ,VZ( 100) ,DT( 100) ,NX( 100) 
DIMENSION  NY( 100) ,NZ( 100) ,RDV( 100) ,EX( 100) ,EY( 100) ,EZ( 100) 
DIMENSION  NDEC 100) ,EDR( 100) ,V( 100) ,HT( 100) ,RDRF( 100) , INTA( 100) 
DIMENSION  R0MA( 100) ,THJO( 100) ,ATE( 100) ,CTE( 100) 

COMMON/OBLATE 1/DAY , RX , RY , RZ , VX , VY , VZ , DT , HX , HY , HZ , NX , NY , NZ , K , KK 
C0MM0N/0BLATE2/RDV,R,V,EX,EY,EZ,MU,NDE,EDR,H,N,E,P, I>OM,W,F 
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C0MM0N/0BLATE3/PI ,EC,M,A,HT,ER,TH,THD,RTD,MD,WD,OMD, IDjECD 
COMMON / 0BLATE4/FD , LINE , J, THF , THFD , IF , IFD , OMF , OMFD , RF , INT, ROM 
COMMON/OBLATE5/RJ , DR , DID , DTHD , DOMD , ESTERR , ACTERR , TERROR , JVER 
COMMON/OBLATE6/RFX , RFY , RFZ , ARC , ARCD , RDRF , DRV , RJ2 , JN, JORBIT 
COMMON / OBLATE  7 / INTA , ROMA , TH JO , ATE , CTE 

CALCULATE  R,V,E,H,N,P,I,OM,W,F,EC,M,HT,TH 

DO  10  J  =  1,  KK 

CALL  CROSS(RX( J) ,RY( J) ,RZ( J) ,VX(J) ,VY( J) ,VZ( J) ,HX( J) ,HY( J) , 
+  HZ(  J) ) 

CALL  CR0SS( 0. ODO , 0. ODO , 1.  ODO , HX( J) ,HY( J) ,HZ(J) , NX( J) , NY( J) , 
+  NZ( J) ) 

CALL  DOT(RX( J) ,RY( J) ,RZ(J) ,VX( J) ,VY( J) ,VZ( J) ,RDV( J) ) 

R( J)=DSQRT( RX( J)*RX( J) +RY( J)*RY( J)+RZ( J)*RZ( J) ) 
V(J)=DSQRT(VX(J)*VX(J)+VY(J)*VY(J)+VZ(J)*VZ(J)) 

EX( J)=((V( J)*V( J)-MU/R( J) )*RX( J)-RDV( J)*VX( J))/MU 
EY( J)=( ( V( J)*V( J) -MU/R( J) )*RY( J) -RDV( J)*VY( J) ) /MU 
EZ(J)=((V(J)*V(J)-MU/R(J))*RZ(J)-RDV(J)*VZ(J))/MU 
CALL  DOT(NX( J) ,NY(J) ,NZ( J) ,EX( J) ,EY( J) ,EZ( J) ,NDE( J)) 

CALL  DOT(EX( J) ,EY( J) ,EZ( J) ,RX(J) ,RY( J) ,RZ( J) ,EDR( J)) 
H(J)=DSQRT(HX( J)*HX( J)+HY( J)*HY( J)+HZ( J)*HZ( J)) 
N(J)=DSQRT(NX(J)*NX(J)+NY(J)*NY(J)+NZ(J)*NZ(J)) 
E(J)=DSQRT(EX(J)*EX(J)+EY(J)*EY(J)+EZ(J)*EZ(J)) 
P(J)=H(J)*H(J)/MU 
I( J)=DACOS(HZ( J)/H( J)) 

OM(J)=DACOS(NX(J)/N( J)) 

W(J)=DACOS(NDE(J)/(N(J)*E(J))) 

F( J)=DACOS(EDR( J)/(E( J)*R( J))) 

IF(NY( J).  LE.  0.  ODO)THEN 
0M( J)=2. ODO*PI-OM( J) 

ENDIF 

IF(EZ( J).  LE.  0. ODO)THEN 
W( J)=2.  ODO*PI-W(J) 

ENDIF 

IF(RDV( J).  LE.  0.  ODO)THEN 
F( J)=2.  ODO*PI-F( J) 

ENDIF 

EC( J)=DACOS( (EC J)+DCOS(F( J) ) )/( 1.  ODO+E( J)*DCOS(F( J) ) ) ) 

IF(F( J).  GE.  PI)THEN 
EC( J)=2.  ODO*PI-EC(J) 

ENDIF 

M(J)=EC(J)-E(J)*DSIN(EC(J)) 

A( J)=P( J)/( 1-E( J)*E( J)) 

AT=(MU/(N( 1)*N( 1) ) )**( i. 0D0/3. ODO) 

HT( J)=R( J)-ER 

RJ=3.  ODO*RJ2*ER"'ER/(  2.  ODO*P(  1)*P(  1)) 

TH(J)=F(J)+W(J) 

THD( J)=TH( J)*RTD 
20  IF(THD( J) . GT. 360. ODO)THEN 

THD( J)=THD( J)-360. ODO 
GOTO  20 
END  IF 

TH(J)=THD(J)/RTD 
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c 
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c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


c 


c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 


CONVERT  ORBITAL  ELEMENTS  TO  DEGREES 


MD( J)=M( J)*RTD 
WD(J)=W(J)*RTD 
OMD ( J ) =0M( J) *RTD 
ID(J)=I(J)*RTD 
ECD(J)=EC(J)*RTD 
FD(J)=F(J)*RTD 
THD(J)^TH(J)*RTD 
CONTINUE 
RETURN 
END 


SUBROUTINE  CROSS 


rc 

•fc 


* 

* 

*  * 


SUBROUTINE  CROSS ( AX , AY , AZ , BX , B Y , BZ , CX , CY , CZ ) 
'IMPLICIT  DOUBLE  PRECISION  (A-I,M-Z) 

CALCULATE  THE  CROSS  PRODUCT  OF  TWO  VECTORS  A  AND  B 

CX=AY*BZ-AZ*BY 

CY=AZ*BX-AX*BZ 

CZ=AX*BY-AY*BX 

RETURN 

END 

*  * 

*  SUBROUTINE  DOT  * 

*  * 

SUBROUTINE  DOT( AX , AY , AZ , BX , BY , BZ , ADB ) 

IMPLICIT  DOUBLE  PRECISION  (A-I,M-Z) 

CALCULATE  THE  DOT  PRODUCT  OF  TWO  VECTORS  A  AND  B 

ADB=AX*BX+AY*BY+AZ*B Z 

RETURN 

END 


SUBROUTINE  PINITIAL 


SUBROUTINE  PINITIAL 

IMPLICIT  DOUBLE  PRECISION  (A-I.M-Z) 

CHARACTER*20  LINE 


HO 


DIMENSION  M( 100) ,MD( 100) ,E( 100) ,W( 100) ,WD( 100) ,0M( 100) ,OMD( 100) 
DIMENSION  I( 100) , ID( 100) ,F( 100) ,FD( 100) ,EC( 100) ,ECD( 100) ,A( 100) 
DIMENSION  R( 100) ,H( 100) ,N( 100) ,TH( 100) ,THD( 100) 

DIMENSION  RF( 100) , IF( 100) , IFD( 100) ,OMF( 100) ,OMFD( 100) ,THF( 100) 
DIMENSION  THFD( 100) ,P( 100) , JORBIT( 100) ,DR( 100) ,DID( 100) ,DTHD( 100) 
DIMENSION  DOMD( 100) ,RX( 100) ,RY( 100) ,RZ( 100) ,RFX( 100) ,RFY( 100) 
DIMENSION  RFZ( 100) ,DRV( 100) ,ARC( 100) ,ARCD( 100) 3 DAY ( 100) ,HX( 100) 
DIMENSION  HY( 100) ,HZ( 100) , VX( 100) , VY( 100) ,VZ( 100) ,DT( 100) ,NX( 100) 
DIMENSION  NY( 100) ,NZ( 100) ,RDV( 100) ,EX( 100) ,EY( 100) }EZ( 100) 
DIMENSION  NDE( 100) ,EDR( 100) ,V( 100) ,HT( 100) ,RDRF( 100) , INTA( 100) 
DIMENSION  ROMA(IOO) ,THJO( 100) ,ATE( 100) ,CTE( 100) 

COMMON/OBLATE 1 /DAY , RX , RY , RZ , VX  3  VY , VZ , DT , HX , HY , HZ , NX  3  NY , NZ , K , KK 
COMMON/ 0BLATE2 /RDV , R , V  5  EX , E Y , E  Z , MU , NDE , EDR , H , N , E , P , I , OM , W , F 
COMMON/ OBLATE  3/ PI , EC , M , A , HT , ER , TH , THD , RTD , MD , WD , OMD , I D , ECD 
COMMON/OBLATE4/FD , LINE , J , THF , THFD , IF , IFD , OMF , OMFD , RF , INT , ROM 
COMMON/ OBLATE5 /RJ , DR , DID , DTHD , DOMD , ESTERR , ACTERR , TERROR , JVER 
COMMON/OBLATE6/RFX , RFY , RFZ  s ARC , ARCD , RDRF , DRV , RJ2 , JN , JORBIT 
COMMON/ OBLATE  7 / INTA , ROMA ,THJ0, ATE , CTE 
C 

C  PRINT  INITIAL  ORBITAL  ELEMENTS 

C 

WRITE(6, '(/)') 

WRITE(2, '(/)') 

WRITE(6,6000)  '  ORBITAL  ELEMENTS' 

WRITE(2,6000)  '  ORBITAL  ELEMENTS' 


WRITEC6,6100) 

LINE 

WRITEC2,6100) 

LINE 

WRITE (6, 6200) 

'M  = 

'  ,MD( 1) j 'N  =  ' >N( 1) , ' E  =  ’ 

,ECl), 

+ 

'W  = 

',WDC1),'0M  =  '  0MD(1),'I  = 

' ,ID(1) 

+ 

'EC  = 

'  ,ECDC 1) , ' A  =  ' jACl), 'R  = 

' ,R( 1) , 

+ 

'HT  = 

'  ,HT( 1) , 'H  =  ',HC1),'F  =  ' 

>FD( 1) , 

+ 

'TH  = 

'  THDd) 

WRITEC2,6200) 

*M  = 

',MD(1),'N  =  ',N(1)5'E  =  ' 

,E(1), 

+ 

•w  = 

'  ,WDC1),'0M  =  ' ,OMD(l) ,  I  = 

MD(1) 

+ 

'EC  = 

'  ,ECD(1).'A  =  ' jA(l) , 'R  = 

',R(D, 

+ 

'HT  = 

',HTC1),'H  =  ' ,H( 1) , 1 F  =  ' 

,FD(1), 

+ 

'TH  = 

' ,THDC 1) 

C 

6000  FORMAT(3X,A) 

6100  FORMATC  3X , A20/ ) 

6200  FORMATC 13(3X,A5,D18. 10/)/) 
C 

RETURN 

END 


C 

C 

C 

C 

C 

C 

C 


*****Vt**A5V***Vf**********Vf****')Wf** 
*  * 


*  SUBROUTINE  INTEGRATE  * 

*  * 


*****>V**5V*5V>VVf*5'c********-V-V******Vf* 


SUBROUTINE  INTEGRATE 

IMPLICIT  DOUBLE  PRECISION  CA-I,M-Z) 

CHARACTER* 20  LINE 

DIMENSION  M( 100) ,MD( 100) ,E( 100) ,W( 100) ,WDC 100) ,0M( 100) ,OMD( 100) 
DIMENSION  I( 100) , ID( 100) ,FC 100) ,FD( 100) ,EC( 100) ,ECD( 100) ,A( 100) 
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DIMENSION  R( 100) ,H( 100) ,N( 100) ,TH( 100) ,THD( 100) 

DIMENSION  RF( 100) , IF( 100) , IFD( 100) ,OMF( 100) ,OMFD( 100) ,THF( 100) 
DIMENSION  THFDC1 00) ,P( 100) , JORBIT( 100) ,DR( 100) ,DID( 100) ,DTHD( 100) 
DIMENSION  DOME  ',0)  ,RX(  100)  ,RY(  100)  ,RZ(  100)  ,RFX(  100)  ,RFY(100) 
DIMENSION  RFZ( J  •-)  ,DRV(  100)  ,ARC(  100)  ,ARCD(  100)  ,DAY(  100)  ,HX(  100) 
DIMENSION  HY( lOo) ,HZ( 100) ,VX( 100) ,VY( 100) ,VZ( 100) ,DT( 100) ,NX( 100) 
DIMENSION  NY( 100) ,NZ( 100) ,RDV( 100) ,EX( 100) ,EY( 100) ,EZ( 100) 
DIMENSION  NDE( 100) ,EDR( 100) ,V( 100) ,HT( 100) ,RDRF( 100) , INTA( 100) 
DIMENSION  ROMA( 100) ,THJO( 100) ,ATE( 100) ,CTE( 100) 
COMMON/OBLATEl/DAY,RX,RYjRZ,^X,VY,VZ jDTjHX,HY,HZ,NX,NY,NZ ,K,KK 
C0MM0N/0BLATE2/RDV,R,V,EXjEY,EZ,MU,NDE,EDR,HjN,E,P, I,OM,W,F 
COMMON/ 0BLATE3/PI , EC , M , A , HT, ER , TH , THD , RTD , MD , WD , OMD , ID , ECD 
C0MM0N/0BLATE4/FD , LINE , J ,THF ,THFD , IF , IFD , OMF , OMFD , RF , INT , ROM 
COMMON/OBLATE5 /RJ , DR , DID , DTHD , DOMD , ESTERR , ACTERR , TERROR , JVER 
COMMON/OBLATE6/RFX , RFY , RFZ , ARC , ARCD , RDRF , DRV , RJ2 , JN , JORB IT 
COMMON/ OBLATE  7 / INTA , ROMA , TH JO , ATE , CTE 

EQUATE  INITIAL  VALUES 

THF( 1)=TH( 1) 

THFD( 1)=THD(1) 

IF( 1)=I( 1) 

IFDC 1)=ID( 1) 

OMFC l)=OM( 1) 

OMFDC 1)=0MD( 1) 

RF( 1)=RC 1) 

ESTIMATE  UPPER  BOUND  OF  THETA 
THF( J)=TH( J)+( J-1)*0. 57D0*2. 0D0*PI 
INITIALIZE 
C 1X10)-12 

ESTERR=0. 000000000001D0 
INT=DT( J)*H( 1) 

DTHF=0.  1745329251994 
C  C1X10)-10 

TERR0R=0.  OOOOOOOOOIDO 
C 

10  CALL  ROMBERG 

ACTERR=INT-ROM 
IF( ACTERR.  LT.  0.  0D0)THEN 
THF ( J ) =THF ( J) -DTHF 
GOTO  10 
ENDIF 

TEMPTHF=THFC  J) 

GOTO  20 
C 

30  CALL  ROMBERG 

ACTERR=INT-ROM 

20  IF(ACTERR.  GE. 0. 0D0)THEN 

IF(C ACTERR/ INT). LE. ESTERR)THEN 
GOTO  40 
ELSE 

TEMPTHF=THF( J) 


ROM00080 


ROM00190 
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THF( J)=THF( J)+DTHF 
GOTO  30 
ENDIF 
ELSE 

DTHF=DTHF/2. ODO 
THF( J)=TEMPTHF+DTHF 
GOTO  30 
END  IF 
C 

40  INTA(J)=INT 
R0MA( J)=R0M 
RETURN 
END 

********************************* 

*  * 

*  SUBROUTINE  ROMBERG  * 

*  * 

********************************* 

SUBROUTINE  ROMBERG  ROM00190 

IMPLICIT  DOUBLE  PRECISION  (A-I.M-Z) 

CHARACTER*20  LINE 

DIMENSION  M( 100) ,MD( 100) ,E( 100) ,W( 100) ,WD( 100) ,OM( 100) ,0MD( 100) 

DIMENSION  I( 100) , ID( 100) ,F( 100) ,FD( 100) ,EC( 100) ,ECD( 100) ,A( 100) 

DIMENSION  R( 100) ,H( 100) ,N( 100) ,TKC 100) ,THD( 100) 

DIMENSION  RF( 100) , IF( 100) , IFD( 100) ,0MF( 100) ,OMFD( 100) ,THF( 100) 

DIMENSION  THFD( 100) ,P( 100) , JORBITC 100) ,DR( 100) ,DID( 100) ,DTHD( 100) 
DIMENSION  D0MD( 100) ,RX( 100) ,RY( 100) ,RZ( 100) ,RFX( 100) ,RFY( 100) 

DIMENSION  RFZ( 100) ,DRV( 100) ,ARC( 100) ,ARCD( 100) ,DAY( 100) ,HX( 100) 

DIMENSION  HY( 100) ,HZ( 100) , VX( 100) ,VY( 100) , VZ( 100) ,DT( 100) ,NX( 100) 
DIMENSION  NY( 100) ,NZC 100) ,RDV( 100) ,EX( 100) ,EYC 100) ,EZ( 100) 

DIMENSION  NDEC 100) ,EDRC 100) , VC  100) ,HTC 100) ,RDRFC 100) , INTAC 100) 

DIMENSION  ROMAC 100) ,THJOC 100) ,ATEC 100) ,CTEC 100) 

COMMON/OBLATE 1/DAY , RX , RY , RZ , VX , VY , VZ , DT , HX , HY , HZ , NX , NY , NZ , K , KK 
COMMON/OBLATE2/RDV ,  R ,  V ,  EX ,  E  Y.,  EZ ,  MU ,  NDE ,  EDR ,  H ,  N ,  E ,  P ,  I ,  OM ,  W ,  F 
COMMON/OBLATE3/PI , EC , M , A , HT , ER , TH , THD , RTD , MD , WD , OMD , ID , ECD 
C0MM0N/0BLATE4/FD , LINE , J , THF , THFD , IF , IFD , OMF , OMFD , RF , INT , ROM 
C0MM0N/0BLATE5 /R J , DR , DID , DTHD , DOMD , ESTERR , ACTERR , TERROR , JVER 
COMMON/OBLATE6/RFX , RFY , RFZ , ARC , ARCD , RDRF , DRV , RJ2 , JN , JORB IT 
COMMON/ OBLATE  7 / INTA , ROMA , TH JO , ATE , CTE 


EXTERNAL  FUNC  R0M00460 

INITIALIZE  VARIABLES  R.OM00560 

HS=THFC  J) -THFC 1 )  ROM00570 

FUNCA=FUNCCRJ, AC  1) , IC 1) ,EC 1) ,WC 1) ,THC 1) ,TKF( 1) , JVER) 
FUNCB=FUNC(RJjAC1),I(1),EC1),WC1),THC1),THFCJ),JVER) 

P  C 1 ) =HS*( FUNCA+FUNCB ) / 2 . ODO  R0M00580 

JM=1  R0M00590 

ROM00610 

BEGIN  THE  ROMBERG  LOOP.  ROM00620 

DO  10  JN  =  1,  100  ROM00630 

OLDINT=P( 1)  ROM00640 
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HH=HS 
SS=0.  0D0 

TT=THF( l)+HH/2. 0D0 
DO  20  L  =  1,  JM 

t-tt 

SS=SS+FUNC ( R J , A( 1) ,  I( 1) ,E( 1) ,W( 1) ,TH( 1) ,T, JVER) 
TT=TT+HH 

20  CONTINUE 

SUM=HH*SS 

P( JN+1)=(P( JN)+SUM)/2. 0D0 
D=l. 0D0 

DO  30  JK  =  JN,  1,  -1 
D=4. 0D0*D 

P( JK)=P( JK+1)+(P( JK+1) -P( JK))/(D-1. 0D0) 

30  CONTINUE 

ERROR=(P(1)-OLDINT) 

IF( JN.  GE.  10)THEN 

IF  (DABS(ERROR/OLDINT).  LE. TERROR) THEN 
GOTO  40 
ENDIF 
ENDIF 

HS=HS/2.  0D0 
JM=JM*2 

10  CONTINUE 

40  ROM=P( 1); 

RETURN 

END 


R0M00670 

R0M00680 

R0M00690 

R0M00720 

ROM00730 

R0M00740 

R0M00750 

R0M00760 

ROM00780 

R0M00800 

ROM00830 

R0M00840 

R0M00850 

R0M00860 

R0M00870 

R0M00900 

ROM00930 


R0M00940 

R0M00950 

R0M00960 

ROM00970 


*A***A***?V****?V;iV***5V********A5V5V5V* 
*  * 


*  FUNCTION  FUNC  * 

*  * 
************* ******************** 


FUNCTI ON  FUNC (RJ,A1,I1,E1,W1,TH1, THFJ , JVER ) 

IMPLICIT  DOUBLE  PRECISION  (A-I,0-Z)  R0M00030 

EXTERNAL  RADIUS 

S=DSIN( II) 

S2=S*S 
S4=S2*S2 
S6=S4*S2 
E2=E1*E1 

RAD=RADIUS(RJ,A1,I1,E1,W1,TH1,THFJ, JVER) 

F  (  SOLUTION  ) 

IF(JVER.  EQ.  1)THEN 
C 

Yl=112.  ODO-75.  0D0*S6+260. 0D0*S4-296. 0D0*S2 
Y 2=R J*THF J* ( 2 .  5D0*S2-2. 0D0) 

Y3=2. 0D0*W1-Y2 

Y4=24.  0D0*(5. 0D0*S2-4. 0D0)*(5. 0D0*S2-4. 0D0) 

Y5=E2*S2*( 14. 0D0-15. 0D0*S2)*(15. 0D0*S2-13.  ODO) 
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Y6=9. 0D0*E2+34. ODO 

Y9=12. 0D0*(5. 0D0*S2-4. ODO) 

Y12=9.  0D0*E2-34. ODO 

YF=(2. 5D0*S2-2.  ODO)*(THFJ-TH1)+E2*Y1*DSIN(Y2)*DCOS(Y3)/Y4 

YS=Y5*DCOS(2.  0D0*Wl)/(2. 0D0*Y9)+ 

+  E1*S2*( 15.  0D0*S2-13.  0D0)*DC0S(THl+Wl)/2.  ODO+ 

+  E1*S2*( 15. 0D0*S2-13.  ODO)*DCOS(3.  0D0*THl-Wl)/6.  ODO+ 

+  S2*( 15.  0D0*S2-13.  0D0)*DC0S(2.  0D0*THl)/2.  ODO+ 

+  (5.  0D0*Y6*S4+4.  ODO*Y12*S2-56.  0D0*E2)/96. ODO 

Y=THF J - W 1+R J*YF+R J*R J*THF J*YS 

F=(2. 0D0-3. 0D0*S2)*DC0S(2.  ODO*THFJ)/2.  ODO+ 

+  El*(S2-l)*DCOS(Y)+El*(3.  ODO-4. ODO*S2)*DCOS(Y+2. 0D0*THFJ)/6. ODO+ 

+  El*( 1. 0D0-2. 0D0*S2)*DC0S(Y-2.  0D0*THFJ)/2.  0D0+S2-1.  ODO+ 

+  E2*S2*( 15. 0D0*S2-14.  0D0)*DSIN(Y2)*DSIN(Y3)/Y9+ 

+  S2*DC0S(2.  ODO*THl) /2.  0D0+E1*S2*DC0S(3. 0D0*THl-Wl)/6. ODO+ 

+  El*S2*DCOS(THl+Wl)/2.  ODO 

ENDIF 

F  (  SIMPLIFIED  SOLUTION  ) 

IF( JVER. EQ.  2) THEN 

F=(2.  0D0-3. ODO*S2)*DCOS(2.  0D0*THFJ)/2.  ODO+ 

+  El*(S2-l)*DCOS(THFJ-Wl)+ 

+  El*(3. ODO-4. 0D0*S2)*DC0S(3.  0D0*THFJ-Wl)/6.  ODO+ 

+  El*(l. 0D0-2. ODO*S2)*DCOS(THFJ+Wl)/2.  0D0+S2-1. ODO+ 

+  S2*DC0S(2.  0D0*THl)/2.  ODO+E1*S2*DCOS(3. 0D0*THl-Wl)/6.  ODO+ 

+  E 1*S  2*DCOS ( TH 1+W1 ) / 2 .  ODO 

ENDIF 

F  (  TWO  BODY  SOLUTION  ) 

IF( JVER. EQ.  3)THEN 
F=0. ODO 

ENDIF 


FUNCTION 

FUNC=RAD*RAD*( 1.  ODO+RJ*F) 


RETURN 

END 


**Vf5V****5'f**5V-.'f**5VVt***5'f************ 
*  * 

*  FUNCTION  RADIUS  * 

*  * 
****  *  *  ***  **  ******  ******  *  *  ******** 


FUNCTION  R ADIUS ( R J , A 1 , I 1 , E 1 , W 1 , TH 1 , TF J , JVER ) 


OGG  O  O  O 


IMPLICIT  DOUBLE  PRECISION  (A-I,0-Z) 

CALCULATE  E,  SINE,  AND  COSINE  FUNCTIONS 

S=DSIN(I1) 

S2=S*S 
S4=S2*S2 
S6=S4*S2 
C=DCOS(Il) 

C2=C*C 

SC=DSIN(I1)*DC0S(I1) 

E2=E1*E1 
P0=A1*(1. 0D0-E2) 

RADIUS  BOTTOM  (  SOLUTION  ) 

IF(JVER. EQ.  1)THEN 
C 

Yl=112.  ODO-75.  0D0*S6+260.  0D0*S4-296.  0D0*S2 

Y2=RJ*TFJ*( 2. 5D0*S2-2.  0D0) 

Y3=2.  0D0*W1-Y2 

Y4=24. 0D0*(5. 0D0*S2-4. 0D0)*(5. 0D0*S2-4.  0D0) 

Y5=E2*S2*( 14. 0D0-15. 0D0*S2)*( 15. 0D0*S2-13. 0D0) 

Y6=9.  0D0*E2+34.  0D0 

Y9=12. 0D0*(5.  0D0*S2-4.  ODO) 

C 

Yll=15. 0D0*(2.  0D0+E2)*S4-14. 0D0*(4.  0D0+E2)*S2+24.  ODO 
C 

Y12=9. 0D0*E2-34.  ODO 
C 

YF=(2. 5D0*S2-2.  ODO)*(TFJrTHl)+E2*Yl*DSIN(Y2)*DCOS(Y3)/Y4 
C 

YS=Y5*DCOS( 2.  0D0*Wl)/(2.  0D0*Y9)+ 

+  E1*S2*( 15. 0D0*S2-13. ODO)*DCOS(THl+Wl)/2.  ODO+ 

+  E1*S2*( 15.  0D0*S2-13.  ODO)*DCOS(3. 0D0*THl-Wl)/6.  ODO+ 

+  S2*( 15. 0D0*S2-13. ODO)*DCOS(2. 0D0*THl)/2.  ODO+ 

+  (5. 0D0*Y6*S4+4. 0D0*Y12*S2-56. 0D0*E2)/96.  ODO 

C 

Y=TF  J - W 1+R J*YF +R J*R J*TF J* YS 
C 

RF1=1. ODO-1.  5D0*S2+E2*( 1.  ODO-1.  25D0*S2)- 
+  ((2.  0D0+5.  0D0*E2)*S2-2. 0D0*E2)*DC0S(2.  0D0*TFJ)/12.  ODO+ 

+  E2*(9. ODO*S2-8. 0D0)*DC0S(2. 0D0*Y)/12.  ODO+ 

+  El*( -11.  0D0*S2+6. ODO)*DCOS(Y+2. 0D0*TFJ)/24.  ODO+ 

+  E2*( -3.  0D0*S2+2.  ODO)*DCOS(2. 0D0*Y+2.  0D0*TFJ)/24.  ODO+ 

+  E2*(3. 0D0*S2-2. 0D0)*DC0S(2. 0D0*Y-2.  0D0*TFJ)/8.  0D0+ 

+  E 1*Y 1 1*DS IN( Y2 ) *DSIN( TFJ+W1) / Y9 

C 

RF2=E2*S2*( 15.  0D0*S2-14.  0D0)*DSIN(Y2)*DSIN(Y3)/(0.  5D0*Y9)- 
+  E2*S2*DC0S(Y-THl+3.  0D0*W1)/16.  ODO+ 

+  E2*(3.  0D0*S2-2. ODO)*DCOS(Y-3. 0D0*THl+3.  0D0*Wl)/24.  D0- 

+  E2*S2*DC0S(Y-5.  0D0*TH1+3..0D0*W1)/16.  ODO+ 

+  El*(3. 0D0*S2-2. ODO)*DCOS(Y-2. 0D0*THl+2.  0D0*Wl)/4.  ODO- 

+  3. ODO*E l*S2*DCOS ( Y-4. 0D0*THl+2.  0D0*Wl)/8.  ODO- 

+  E1*(S2+1. ODO)*DCOS(Y+2. 0D0*Wl)/4. ODO 

C 


RQM00030 
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RF3=((5.  0D0*E2-2.  0D0)*S?.-2.  0D0*E2)*DC0S(Y+THl+Wl)/8.  0D0+ 

+  ((5. 0D0*E2+6. 0D0)*S2-4. 0D0*(E2+1. 0D0) )*DC0S(Y-THl+Wl)/4.  0D0+ 

+  (2.  0D0*E2-S2*(5.  0D0*E2+14.  0D0))*DC0S(Y-3.  0D0*THl+Wl)/24- 0D0+ 

+  E2*(9. 0D0*S2-4.  0D0)*DC0S(Y+3. 0D0*THl-Wl)/48.  D0+ 

+  E2*(6. 0D0-7.  ODO*S2)*DCOS(Y+THl-Wl)/8.  0D0+ 

+  E2*(4..0D0-5.  0D0*S2)*DC0S(Y-TH1-W1)/ 16.  0D0+ 

+  El*(2.  0D0*S2-1. 0D0)*DC0S(Y+2. 0D0*THl)/4. 0D0 

RF4=E1*( 1. 0D0-3.  0D0*S2)*DC0S(Y-2.  0D0*THl)/4.  0D0+ 

+  El*(2. ODO-3.  0D0*S2)*DC0S(Y)/4. 0D0+ 

+  E1*S2*DC0S(TH1+W1)+S2*DC0S(2. 0D0*TH1)+ 

+  E1*S2*DC0S(3.  0D0*THl-Wl)/3.  0D0 

RFB=1. ODO+E1*DCOS(Y)+RJ*(RF1+RF2+RF3+RF4) 

ENDIF 

RADIUS  BOTTOM  (  SIMPLIFIED  SOLUTION  ) 

IFCJVER.  EQ.  2)THEN 

RF1=1.  0D0-1.  5D0*S2+E2*( 1.  0D0-1.  25D0*S2)- 
+  ((2.  0D0+5. 0D0*E2)*S2-2.  0D0*E2)*DC0S(2.  0D0*TFJ)/12.  0D0+ 

+  E2*(9.  0D0*S2-8.  0D0)*DC0S(2. 0D0*(TFJ-W1) )/12. 0D0+ 

+  El*(-ll.  0D0*S2+6.  ODO)*DCOS(3.  0D0*TFJ-Wl)/24.  0D0+ 

+  E2*( -3.  0D0*S2+2.  ODO)*DCOS(4. 0D0*TFJ-2. 0D0*tfl)/24. 0D0+ 

+  E2*(3. 0D0*S2-2. 0D0)*DC0S(2. ODO*Wl)/8.  0D0- 

+  E2*S2*DC0S(TFJ-THl+2.  0D0*W1)/16.  0D0 

RF2=E2*(3. 0D0*S2-2.  ODO)*DCOS(TFJ-3. 0D0*THl+2. 0D0*Wl)/24.  D0- 
+  E2*S2*DCOS(TFJ-5.  0D0*THl+2.  0D0*W1)/16.  0D0+ 

+  El*(3. 0D0*S2-2.  0D0)*DC0S(TFJ-2.  0D0*THl+Wl)/4.  0D0- 

+  3. ODO*El*S2*DCOS(TFJ-4.  0D0*THl+Wl)/8.  0D0- 

+  E1*(S2+1. 0D0)*DC0S(TFJ+Wl)/4.  0D0+ 

+  ((5. ODO*E2-2. 0D0)*S2-2. 0D0*E2)*DC0S(TFJ+THl)/8.  ODO+ 

+  ((5. ODO*E2+6. 0D0)*S2-4. 0D0*(E2+1. 0D0) )*DCOS(TFJ-THl)/4.  DO 

RF3=( 2. 0D0*E2-S2*(5. 0D0*E2+14. 0D0) )*DC0S(TFJ-3. D0*THl)/24.  D0+ 

+  E2*(9.  0D0*S2-4.  0D0)*DC0S(TFJ+3. 0D0*THl-2. 0D0*Wl)/48. D0+ 

+  E2*(6.  0D0-7.  0D0 *S  2 ) *DCOS ( TF J+TH 1 - 2 .  0D0*Wl)/8.  0D0+ 

+  E2*(4. 0D0-5.  0D0*S2)*DC0S(TFJ-THl-2.  0D0*W1)/16.  0D0+ 

+  El*(2.  0D0*S2-1.  ODO)*DCOS(TFJ+2. 0D0*THl-Wl)/4. 0D0+ 

+  El*( 1.  ODO-3.  ODO*S2)*DCOS(TFJ-2.  ODO*THl-Wl)/4. 0D0+ 

+  El*(2.  ODO-3.  ODO*S2)*DCOS(TFJ-Wl)/4.  0D0 

RF4=E 1*S  2*DCOS ( TH 1+W 1 ) +S  2*DC  0  S ( 2 . OD  0*TH 1 ) + 

+  El*S2*DCOS(3. 0D0*THl-Wl)/3.  0D0 

RFB=E1*DC0S(TFJ-W1+RJ*(TFJ-TH1)*(2.  5D0*S2-2.  0D0))+ 

+  1. ODO+R J*( RF 1+RF2+RF3+RF 4 ) 


ENDIF 

RADIUS  BOTTOM  (  TOO  BODY  SOLUTION  ) 
IF(JVER. EQ.  3)THEN 
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RFB=1.  ODO+E1*DCOS(TFJ-W1) 

END  IF 

RADIUS 

RADIUS=PO/RFB 

RETURN  ROM00410 

END  ROM00420 

******************************* 

*  * 

*  SUBROUTINE  FORMULA  * 

*  * 

******************************* 


SUBROUTINE  FORMULA 

IMPLICIT  DOUBLE  PRECISION  (A-I,M-Z) 

CHARACTER*20  LINE 

DIMENSION  M( 100) ,MD( 100) ,E( 100) ,W( 100) ,WD( 100) ,0M( 100) ,0MD( 100) 
DIMENSION  I( 100) ,  ID( 100) ,F( 100) ,FD( 1 00) ,EC( 100) ,ECD( 100) , A( 100) 
DIMENSION  R( 100) ,H( 100) ,N( 100) ,TH( 100) ,THD( 100) 

DIMENSION  RF( 100) j IF( 100) , IFD( 100) ,0MF( 100) ,0MFD( 100) ,THF( 100) 
DIMENSION  THFD(  100)  ,P(  l00) , JORBIT(  100)  ,DR(  100)  ,DID(  100)  ,DTHD(  100) 
DIMENSION  D0MD( 100) ,RX( 100) ,RY( 100) ,RZ( 100) ,RFX( 100) ,RFY( 100) 
DIMENSION  RFZ( 100) ,DRV( 100) ,ARC( 100) ,ARCD( 100) ,DAY( 100) ,HX( 100) 
DIMENSION  HY( 100) ,HZ( 100) ,VX( 100) ,VY( 100) ,VZ( 100) ,DT( 100) ,NX( 100) 
DIMENSION  NY( 100) ,NZ( 100) ,RDV( 100) ,EX( 100) ,EY( 100) ,EZ( 100) 
DIMENSION  NDE( 100) ,EDR( 100) , V( 100) ,HT( 100) ,RDRF( 100) , INTA( 100) 
DIMENSION  ROMAC 100) ,THJO( 100) ,ATE( 100) ,CTE( 100) 

COMMON/OBLATE 1/DAY , RX , RY , RZ , VX , VY , VZ , DT , HX , HY , HZ , NX , NY , NZ , K , KK 
COMMON/ 0BLATE2/RDV , R , V , EX , EY ,EZ , MU , NDE , EDR , H , N , E , P , I , OM , W , F 
C0MM0N/0BLATE3/PI , EC , M , A , HT , ER , TH , THD ,RTD } MD , WD , OMD , ID , ECD 
C0MM0N/0BLATE4/FD , LINE , J , THF , THFD , IF , IFD , OMF , OMFD , RF , INT , ROM 
COMMON/OBLATE5 /RJ , DR , DID , DTHD , DOMD , ESTERR , ACTERR , TERROR , JVER 
COMMON/OBLATE6/RFX , RFY , RFZ , ARC , ARCD , RDRF , DRV , RJ2 , JN , JORB IT 
COMMON/OBLATE  7/ INTA , ROMA , THJO , ATE , CTE 

EXTERNAL  RADIUS 

CALCULATE  E,  SINE,  AND  COSINE  FUNCTIONS 

S=DSIN( I( 1) ) 

S2=S*S 

S4=S2*S2 

S6=S4*S2 

C=DCOS(I(l)) 

C2=C*C 

SC=DSIN( I ( 1 ) ) *DC0S ( I ( 1 ) ) 

E2=E( 1)*E( 1) 


FORMULA  (  SOLUTION  ) 

IF( JVER. EQ. 1)THEN 
C 

Y1=U2.  ODO-75.  0D0*S6+260.  0D0*S4-296.  0D0*S2 


1  IS 


ooooooo  oooo  ooo 


Y2=RJ*THF(  J)*(2. 5D0*S2-2.  0D0) 

Y3=2.  0D0*W(1)-Y2 

Y4=24.  0D0*(5.  0D0*S2-4.  0D0)*(S.  0D0*S2-4. 0D0) 

Y5=E2*S2*( 14.  0D0-15.  0D0*S2)*(15.  0D0*S2-13.  0D0) 

Y6=9.  0D0*E2+34.  0D0 

Y7=15.  0D0*S4-45. 0D0*S2+28.  ODO 

Y8=6.  0D0*(5.  ODO*S2-4.  0D0)*(5.  0D0*S2-4.  ODO)  a 

Y9=12.  0D0*(5.  0D0*S2-4.  ODO) 

C 

Y10=(6.  0D0-S2)/12.  0D0-E( 1)*S2*DC0S(3.  0D0*TH(l)-W(l))/3.  0D0- 
+  S2*DC0S(2.  0D0*TH( 1))+E2*(7. 0D0*S2-4. 0D0)/24. ODO 

C 

'Y12=9.  0D0*E2-34.  ODO 
C 

YF=( 2. 5D0*S2-2. 0D0)*(THF( J) -TH( 1) )+E2*Yl*DSIN(Y2)*DC0S(Y3)/Y4 
C 

YS=Y5*DC0S(  2.’0D0*W(  1)  )  /(  2.,ODO*Y9)+ 

+  E( 1)*S2*( 15. 0DQ*S2-13. ODO)*DCOS(TH(l)+W(l))/2.  ODO+ 

+  E( 1)*S2*( 15.  0D0*S2-13.  0D0)*DC0S(3. 0D0*TH( 1) -W( 1) )/6. 0D0+ 

+  S2*( 15.  0D0*S2-13.  ODO)*DCOS(2. 0D0*TH( 1) )/2.  0D0+ 

+  (5.  0D0*Y6*S4+4.  0D0*Y12*S2-56. 0D0*E2)/96. ODO 

C. 

Y=THF( J) -W(1)+RJ*YF+RJ*RJ*THF( J)*YS 


CALCULATE  INCLINATION  (  SOLUTION  ) 

IF(J)=I(1)+SC*RJ*(DC0S(2.  0D0*THF( J))/2. 0D0+ 

+  E(1)*DCOS(Y+2.0DO*THF(J))/6.DO+E(1)*DCOS(Y-2.DO*THF(J))/2.DO+ 

+  E2*(14.  DO-15.  D0*S2)*DSIN(Y2)*DSIN(Y3)/( 12. D0*(5. D0*S2-4. DO))- 

+  DC0S(2. 0D0*TH( l))/2.  0D0-E( 1)*DC0S(3.  0D0*TH( l)-W(l))/6.  0D0- 

+  E( l)*DCOS(TH( 1)+W( 1) )/2.  ODO) 

CALCULATE  LONGITUDE  OF  THE  ASCENDING  NODE  (  SOLUTION  ) 


RJ2=RJ*RJ 


OMF( J)=OM( 1)+C*RJ*(TH( 1) -THF( J)+DSIN(2. 0D0*THF( J) )/2. 0D0- 
+  E(l)*DSIN(Y)+E(l)*DSIN(Y+2.  0D0*THF( J))/6.  ODO- 

+  E(l)*DSIN(Y-2.  0D0*THF( J) )/2.  0D0-DSIN(2.  0D0*TH( 1) )/2.  0D0+ 

+  E( 1)*DSIN(TH( 1)-W( 1) ) -E( 1)*DSIN(3.  0D0*TH( 1)-W( l))/6. 0D0- 

+  E( 1)*DSIN(TH( 1)+W( 1 ) ) / 2.  ODO+E2*Y7*DSIN(Y2)*DCOS(Y3)/Y8)+ 

+  C*R J2*THF( J)*( E2*S2*( 15.  ODO*S2-14.  ODO)*DCOS(2. D0*W( 1))/Y9- 

+  E( l)*S2*DCOS(TH( 1)+W( 1))+Y10) 

END  IF 

FORMULA  (  SIMPLIFIED  SOLUTION  ) 

IF( JVER. EQ. 2) THEN 

CALCULATE  INCLINATION  (  SIMPLIFIED  SOLUTION  ) 

IF(J)=I(l)+SC*RJ*(DCOS(2. ODO*THF( J) )/2. 0D0+ 

+  E(l)*DCOS(3. 0D0*THF(J)-W(l))/6.  0D0+ 

+  E(l)*DC0S(THF(J)+W(l))/2. ODO-DCOS(2.  0D0*TH( 1) )/2. 0D0- 

+  E(l)*DCOS(3. 0D0*TH( 1) -W( 1) )/6. ODO- 
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E(:1)*DC0S(,TH(  1)+W(  l))/2.  ODO) 


C 

C 

C 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 


10 


20 


C 

C 

C 


CALCULATE  LONGITUDE  OF  THE  ASCENDING  NODE  (  SIMPLIFIED  SOLUTION  ) 

OMF(  J)=OM(  1)+C*EJ*(TH( l)'-fHF(  J)+DSIN( 2.  ODO*THF( J) )/2.  ODO- 
+  E(l)*DSINvTHF(J)-W(l))+E(l)*DSIN(3.  ODO*THF( J) -W( 1) )/6. 0D0+ 

^  E(l)*DSIN(THF(J)+W(l))/2. 0D0-DSIN(2. ODO*TH( 1) )/2.  0D0+ 

+  E( 1)*DSIN(TH( 1) -W( 1) ) -E( 1)*DSIN( 3.  ODO*TH( 1) -W( 1) ) /6.  ODO- 

+  E(l)*DSIN(TH(l)+W(l))/2.  ODO) 


ENDIF 

FORMULA  (  TWO  BODY  SOLUTION  ) 

IF(  JVER.  EQ.  3-) THEN 

CALCULATE  INCLINATION  (  TWO  BODY  SOLUTION  ) 

IF( J)=I( 1) 

CALCULATE  LONGITUDE  OF  THE  ASCENDING  NODE  (  TWO  BODY  SOLUTION  ) 
GMF( J)=OM( 1) 

ENDIF 

CALCULATE  RADIUS  (  SOLUTION,  SIMPLIFIED,  OR  TWO  BODY  ) 

RF( J)=RADIUS(RJ,A( 1) ,I( 1) ,E( 1) ,W( 1) ,TH( 1) ,THF( J) , JVER) 

CONVERT  ANGLES  TO  DEGREES 

OMFD(J)=OMF(J)*RTD 
IFD(J)=IF(J)*RTD 
THFD ( J ) =THF ( J ) *RTD 
JORBIT( J)=0 

IF(THFD( J).  GT.  360. ODO)THEN 
THFD(J)=THFD(J)-360.  ODO 
JORBIT( J)=JOR3IT( J)+l 
GOTO  10 
END  IF 

THJO( J)=JORBIT( J)*2.  ODO*PI+TH( J) -TH( 1) 

IF(OMFD( J).  GT.  360. 0D0)THEN 
OMFD( J)=OMFD( J)-360.  ODO 
GOTO  20 
END  IF 

THF( J)=THFD( J)/RTD 
OMF ( J ) =OMFD( J ) /RTD 

CALCULATE  DELTAS 

DR(J)=RF( J)-R( J) 

DID( J)=IFD( J)-ID( J) 

DTHDC  J)=THFD( J) -THD( J) 

IF(DABS(DTHD( J)). GE. 180. ODO) THEN 
IF(DTHD( J).  LT,  0. ODO) THEN 
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DTHD( J)=DTHD( J)+360.  0D0 
ELSE 

DTHD( J)=DTHD( J)-360.  0D0 
ENDIF 
ENDIF 

D0MD(  J.)=0MFD(  J)  -OMD(  J) 

RETURN 

END 


* 

* 

* 

SUBROUTINE  INERTIAL 

* 

* 

* 

Vr******VtA5Wr#r********5'oWr**TV***5V* 


SUBROUTINE  INERTIAL 

IMPLICIT  DOUBLE  PRECISION  (A-I,M-Z) 

CHARACTER*20  LINE 

DIMENSION  M( 100) ,MD( 100) ,E( 100) SW( 100) ,WD( 100) ,0M( 100) ,0MD( 100) 
DIMENSION  I(  100)  , ID(  100)  ,F.(  100)  ,FD(  100)  ,EC(  100)  ,ECD(  100)  ,A(  100) 
DIMENSION  R(  100)  ,H(  100)-jN(  100)  ,TH(  100)  ,THD(  100) 

DIMENSION  RF( 100) , IF(100) , IFD( 100) ,0MF( 100) ,OMFD( 100) ,THF(100) 
DIMENSION  THFD( 100) ,P( 100) , JORBITC 100) ,DR( 100) ,DID( 100) ,DTHD( 100) 
DIMENSION  D0MD( 100) ,RX( 100) ,RY( 100) ,RZC 100) ,RFX( 100) ,RFY( 100) 
DIMENSION  RFZ( 100) ,DRV( 100) ,ARC( 100) ,ARCD( 100) ,DAY( 100) ,HX( 100) 
DIMENSION  HY( 100) ,HZ( 100) ,VX( 100) , VY( 100) ,VZ( 100) ,DT( 100) ,NX( 100) 
DIMENSION  NY( 100) ,NZC 100) ,RDVC 100) ,EXC 100) ,EYC 100) ,EZC 100) 
DIMENSION  NDEC 100) ,EDR( 100) , VC  100) ,HT( 100) ,RDRFC 100) , INTAC 100) 
DIMENSION  EARCC 100) ,EARCDC 100) ,PDRC 100) 

DIMENSION  ROMAC 100) ,THJ0C 100) ,ATE( 100) ,CTEC 100) 

COMMON/OBLATE 1/DAY, RX,RY,RZ,VX,,/Y,VZ>DT,HX,HY, HZ, NX, NY, NZ,K,KK 
C0MM0N/0BLATE2/RDV,R,V,EX,EY,EZ,MU,NDE,EDR,H,N,E,P,I,0M,W,F 
COMMON/ OBLATE  3 / P I , EC , M , A , HT, ER ,TH , THD , RTD , MD , WD , OMD , ID , ECD 
C0MM0N/0BLATE4/FD , LINE , J ,THF ,THFD , IF ,TFD, OMF , OMFD , RF , INT , ROM 
C0MM0N/0BLATE5/RJ , DR , DID , DTHD , DOMD , ESTERR , ACTERR , TERROR , JVER 
COMMON/OBLATE6/RFX , RFY , RFZ , ARC , ARCD , RDRF , DRV , RJ2 , JN , JORBIT 
COMMON/OBLATE7/INTA , ROMA ,THJO , ATE , CTE 
COMMON/ SPECI AL/EARC , E ARCD , PDR , ENG , ENGF 

CALCULATE  INITIAL  ENERGY 

ENG=VCl)*VCl)/2.  ODO-MU/RC 1) 

ENGF=VC  1)*VC -l)/2.  ODO-MU/RFC 1) 

CALCULATE  INERTIAL  COORDINATES 

RFXCJ)=RFCJ)*CDCOSCTHFCJ))*DCOSCOMFCJ))- 
+  DSINCTHFCJ))*DCOSClFCJ))*DSINCOMFCJ))) 

RFYC J)=RFC J)*CDCOSCTHFC J) )*DSINCOMFC J) )+ 

+  DSINCTHFC J) )*DCOS( IFC J) )*DCOSCOMFC J) ) ) 

RFZCJ)=RFCJ)*CDSIN(THFCJ))*DSIN(IF(J))) 

CALCULATE  DR  VECTOR 

DRVC J)=DSQRTC  CRFXC J) -RXC J) )*CRFXC J) -RXC J) )+ 
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+  (RFY(J)-RY(J))*(RFY(J)-RY(J))+ 

+  (RFZ(J)-RZ(J))*(RFZ(J)-RZ(J))) 

PDR(J)=DRV(J)/R(J) 

CALCULATE  ANGLE  ERROR 

CALL  D0T(RX( J) ,RY( J) ,RZ( J) ,RFX( J) ,RFY( J) ,RFZ( J) ,RDRF( J) ) 
ARC(J)=DACOS(RDRF(J)/(R(J)*RF(J))) 

CC=RF( J) 

CCP=R( J) 

BB=ER 

AA=DSQRT(BB*BB+CC*CC-2. ODO*BB*CC*DCOSCARC( J)/2. 0D0)) 
AAP=DSQRT(BB*BB+CCP*CCP-2.  ODO*BB*CCP*DCOS( ARC( J) /2. 0D0) ) 
CCA=PI-DASIN(CC*DSIN(ARC( J)/2.  0D0)/AA) 
CCPA=PI-DASIN(CCP*DSIN(ARC(J)/2.  0D0)/AAP) 

EARC( J)=2.  0D0*PI-CCA-CCPA 
ARCD( J)=ARC( J)*RTD 
E ARCD ( J ) =E ARC ( J ) *RTD 

CALCULATE  DOWNRANGE  AND  CROSSRANGE  ERRORS 


ATE ( J) =R( J) *( DTHD( J) /RTD+DCOS ( I ( J ) )*D0MD( J) /RTD ) 
CTE ( J )=R( J) *( DS IN( TH( J ) )*DID( J) /RTD - 
+  DCOS(TH( J) )*DSIN(I( J) )*DOMD( J) /RTD) 


RETURN 

END 


*5V*5Wf5'c*5V**5V**5VA*5'r**>Wlr*5V!V)VA****)VjV 
*  * 


*  SUBROUTINE  RESULTS 

* 


«v 
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SUBROUTINE  RESULTS 

IMPLICIT  DOUBLE  PRECISION  (A-I,M-Z) 

CHARACTER*20  LINE 
CHARACTER* 11  VERSION 

DIMENSION  M( 100) ,MD(100) ,E( 100) ,W( 100) ,WD( 100) ,OM( 100) ,0MD( 100) 
DIMENSION  I( 100) , ID( 100) ,F( 100) ,FD( 100) ,EC( 100) ,ECD( 100) , A( 100) 
DIMENSION  R( 100) ,H( 100) ,N( 100) ,TH( 100) ,THD( 100) 

DIMENSION  RF( 100) , IFC 100) , IFDC 100) ,OMFC 100) ,OMFDC 100) ,THF( IOC) 
DIMENSION  THFDC 100) ,PC 100) , JORBITC 100) ,DRC 100) ,DIDC 100) ,DTHD( 100) 
DIMENSION  DOMDC 100) ,RXC 100) ,RYC 100) ,RZC 100) ,RFXC 100) ,RFYC 100) 
DIMENSION  RFZC 100) ,DRV( 100) , ARC( 100) , ARCDC 100) ,DAYC 100) ,HXC 100) 
DIMENSION  HY( 100) ,HZC 100) ,VXC 100) ,VYC 100) , VZC 100) ,DTC 100) ,NXC 100) 
DIMENSION  NYC  100) ,NZC ICO) ,RDVC 100) ,EXC 100) ,EYC 100) ,EZ( 100) 
DIMENSION  NDEC 100) ,EDRC 100) ,VC 100) ,HTC 100) ,RDRFC 100) , INTAC 100) 
DIMENSION  EARCC 100) ,E ARCDC 100) ,PDRC 100) 

DIMENSION  ROMAC 100) ,THJOC 100) , ATEC 100) ,CTEC 100) 

COMMON/OBLATE 1/DAY , RX , RY , RZ , VX , VY , VZ , DT , HX , HY , HZ , NX , NY , NZ , K , KK 
C0HM0N/0BLATE2/RDV,R,V,EX,EY,EZ,MU,NDE,EDR,H,N,E,P,I,0M,W,F 
C0MM0N/0BLATE3/PI , EC , M , A , HT , ER , TH,THD , RTD , MD , WD , OMD , ID , ECD 
C0MM0N/0BLATE4/FD , LINE , J , THF , THFD , IF , IFD , OMF , OMFD , RF , INT , ROM 
C0MM0N/0BLATE5/RJ , DR , DID , DTHD , DOMD , ESTERR , ACTERR , TERROR , JVER 
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C0MM0N/0BLATE6/RFX , RFY , RFZ , ARC , ARCD , RDRF , DRV , RJ2, JN , JORB IT 
COMMON/ 0BLATE7/ INTA , ROMA, TH JO , ATE , CTE 
COMMON/ SPEC I AL/EARC , EARCD , PDR , ENG , ENGF 

DR(1)=0.  0D0 
DID( 1)=0. 0D0 
DTHD( 1)=0.  0D0 
D0MD(1)=0.  0D0 
DRV(1)=0.  0D0 
ARCD( 1)=0.  0D0 
EARCD( 1)=0.  0D0 
PDR( 1)=0.  0D0 
THJO( 1)=0.  ODO 
ATE( 1)=0.  ODO 
CTE( 1)=0.  ODO 

IF(  JVER-.  EQ.  1)THEN 
VERSION=' SOLUTION  * 

ELSEIF( JVER.  EQ.  2)THEN 
VERS ION=' SIMPLIFIED' 

ELSEIF( JVER. EQ. 3) THEN 
VERSION^ SECULAR  ' 

ENDIF 

IF(RJ. EQ. 0. ODO)THEN 
VERS ION=' TWO  BODY  ' 

ENDIF 

OUTPUT  RESULTS  FOR  DISSPLA 

IF( JVER.  EQ.  1)THEN 
WRITE(3,3000)  K 
WRITE(3,3100)  RJ 
ENDIF 

DO  10  J  =  I,  KK 

WRITE(3,3100)  DR( J) ,DID(J) ,DTHD( J) 

WRITE(3,3100)  DOMD( J) ,DRV( J) ,EARCD( J) 

WRITE(3,3100)  PDR( J) ,ATE( J) ,CTE( J) 

WRITE(4,3100)  THJO(J) 

CONTINUE 

PRINT  RESULTS 

WRITE(6, '(/)') 

WRITE(2, '(/)') 

WRITE(6,6000)  '  RESULTS' 

WRITE(2,6000)  '  RESULTS' 

WRITE(6,6100)  LINE 
WRITE(2,6100)  LINE 
WRITE(6,6200)  'J  =  ' ,RJ 
WRITE(2,6200)  'J  =  ' ,RJ 
WRITE(6 ,6300)  'VERSION  =  ' , VERSION 
WRITE(2,6300)  'VERSION  =  ' , VERSION 
WRITE(6,6100)  LINE 
WRITE(2,6100)  LINE 
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DO  20  J  =  K,  KK 
C 

WRITE(6,6400)  ’POINT  =  J, 'ORBIT  =  ’,JORBIT(J), 

+  ’ROMBERG  ITERATIONS  =  ' , JN 

WRITE(2,6400)  ’POINT  =  J, 'ORBIT  =  ’,JORBIT(J), 

+  'ROMBERG  ITERATIONS  =  ’ ,JN 

C 

WRITE(6,6500)  *R  =  ’,R(J),'RF  =  ',RF(J),'DR  =  ’ ,DR( J) 

WRITE(2,6500)  'R  =  ',R(J),’RF  =  ',RF(J),’DR  =  ' ,DR(J) 

C 

WRITE(6,6500)  'I  =  '.IDCJ^’IF  =  ’,IFD(J),’DI  =  ’ ,DID( J) 

WRITE(2,6500)  ’I  =  ’,ID(J),’IF  =  ’,IFD(J),’DI  =  ’ ,DID( J) 

C 

WRITE(6,6500)  ’TH  =  ’ ,THD( J) , ’THF  =  ’ ,THFD(J) , ’DTH  =  ’ ,DTHD( J) 
WRITE(2,6500)  ’TH  =  ’ ,THD( J) , ’THF  =  ’ ,THFD( J) , ’DTH  =  ’ ,DTHD(J) 
C 

VrRITE(6,6500)  ’OM  =  ’ ,OMD(J) , *OMF  =  ’ ,OMFD( J) , ’DOM  =  ’ .DOMDCJ) 
WRITE(2,6500)  ’OM  =  ’ ,OMD(J) , ’OMF  =  ’ ,OMFD( J) , 'DOM  =  ’ ,DOMD(J) 
C 

WRITE(6,6500)  ’RX  =  ’,RX(J),’RY  =  ’,RY(J),’RZ  =  ’ ,RZ(J) 
WRITE(2,6500)  ’RX  =  ’,RX(J),'RY  =  ',RY(J),'RZ  =  ’ ,RZ( J) 

C 

WRITE(6,6500)  'RFX  =  ' ,RFX( J) , 'RFY  =  ' ,RFY( J) , 'RFZ  =  ' ,RFZ( J) 
WRITE(2,6500)  'RFX  =  ' ,RFX( J) , ’RFY  =  ' ,RFY( J) , 'RFZ  =  ’ ,RFZ( J) 

C 

WRITE(6,6500)  'DRV  =  ' ,DRV(J) , 'PDR  =  ’,PDR(J), 

+  ’EARC=  ' ,EARCDCJ) 

WRITE(2,6500)  'DRV  =  ' ,DRV(J) , ’PDR  =  ’,PDR(J), 

+  ' EARC=  ' ,EARCD(J) 

C 

WRITE(6,6500)  ’RTE  =  ',DR(J),'ATE  =  ',ATE(J), 

+  'CTE  =  ’ ,CTE(J) 

WRITE(2,6500)  ’RTE  =  ’,DR(J),'ATE  =  ’,ATE(J), 

+  'CTE  =  ' ,CTE(J) 

C 

WRITE(6,6600)  'INT  =  ' ,INTA(J) , ’ROM  =  ' ,ROMA( J) 

KRITE(2,6600)  'INT=  ' ,INTA( J) , 'ROM  =  ’ ,ROMA( J) 

C 

20  CONTINUE 
C 

WRITE(6,6500)  'EG  =  ',ENG,’EGF  =  ' ,ENGF 
WRIlt:C2,6500)  ’EG  =  '  ,ENG,  ’EGF  =  '  ,ENGF 
C 

WRITE(6, ’ (/) ' ) 

WRITE(2, '(/)') 

C 

3000  FORMAT( 3X,I3) 

3100  F0RMAT(3(3X,D18. 10)) 

C 

6000  F0RMAT(3X,A) 

6100  FORMAT( 3X , A20/ / ) 

6200  F0RMAT(3X,A,F8. 6) 

6300  FORMAT( 3X,A,A11//) 

6400  FORMAT(2(3X,A8,I3/) ,3X,A21,I3//) 

6500  FORMAT(3(3X,A6,F23. 15/)) 
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6600  F0RMAT(3(3X,A6,F23. 8/)) 

C 


RETURN 

END 
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